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ABSTRACT 

We present a radiative transfer model, which is apphcable to the study of submil- 
Hmetre spectral line observations of protostellar envelopes. The model uses an exact, 
non-LTE, spherically symmetric radiative transfer 'Stenholm' method, which numer- 
ically solves the radiative transfer problem by the process of 'A-iteration'. We also 
present submillimetre spectral line data of the Class protostars NGC1333-IRAS2 
and Serpens SMM4. We model the data using the Stenholm code. 

We examine the physical constraints which can be used to limit the number and 
range of parameters used in protostellar envelope models, and identify the turbulent 
velocity and tracer molecule abundance as the principle sources of uncertainty in the 
radiative transfer modelling. We explore the trends in the appearance of the predicted 
line profiles as key parameters in the models are varied, such as infall velocity pro- 
file, turbulence and rotation. The formation of the characteristic asymmetric double- 
peaked line profile in infalling envelopes is discussed. We find that the separation of 
the two peaks of a typical infall profile is dependent not on the evolutionary status of 
the collapsing protostar, but on the turbulent velocity dispersion in the envelope. We 
also find that the line shapes can be significantly altered by rotation. 

Fits are found for the observed line profiles of IRAS2 and SMM4 using plausi- 
ble infall model parameters. The density and velocity profiles in our best fit models 
are inconsistent with a singular isothermal sphere model (SIS), since for both objects 
modelled, the infall velocities appear further advanced than a SIS model would pre- 
dict, given the density profile. We find better agreement with a form of collapse which 
assumes non-static initial conditions in agreement with other recent findings. We also 
find some evidence that the infall velocities are retarded from free-fall towards the cen- 
tre of the cloud, probably by rotation, and that the envelope of SMM4 is rotationally 
flattened. 
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1 INTRODUCTION 



The process of star formation is currently not well under- 
stood. However, the main protostellar collapse phases of low- 
mass stars (~O.5-2M0) have been identified observationally 
and labelled as Class (Andre, Ward- Thompson & Barsony 
1993), and Class I (Lada & Wilking 1984; Lada 1987) pro- 
tostars. These are believed to represent the phases during 
which the circumstellar envelope accretes onto the central 
protostar and disk (e.g. Andre 1994; Ward-Thompson 1996). 
The final pre-main-sequence stages of Classes II & III (Lada 
& Wilking 1984; Lada 1987) correspond to the Classical T 
Tauri (CTT) and Weak-line T Tauri (WTT) stages respec- 
tively (Andre & Montmerle 1994). These protostellar stages 



are beginning to be understood at least in outline (for a 
review, see: Andre, Ward- Thompson & Barsony 2000). 

Protostellar infall has been reported in Class sources 
by a number of authors (e.g. Zhou et al. 1993; Ward- 
Thompson et al. 1996) . However, the manner of the collapse 
remains a matter for debate. The ideas of static initial con- 
ditions for collapse (e.g. Shu 1977) have been disputed by 
many authors (e.g. Foster & Chevalier 1993; Whitworth et 
al. 1996). There is now a growing body of evidence which 
suggests that the collapse occurs from non-static initial con- 
ditions, and at a non-constant accretion rate which decreases 
with time (e.g. Kenyon & Hartmann 1995; Henriksen, Andre 
& Bontemps 1997; Safier, McKee & Stabler 1997; Whitworth 
& Ward- Thompson 2001). Infall has even been detected in 
some starless cores, such as L1544 (Tafalla et al. 1999). 
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The densities and temperatures in the gas envelopes 
surrounding the youngest protostars are favourable for ex- 
citing a number of rotational molecular transitions, observ- 
able in the submillimetre waveband. The line profiles of 
these transitions contain information about both the physi- 
cal state and dynamics of the envelope gas, and may poten- 
tially be used to test theoretical models of star formation, 
as many workers have shown (e.g. Bernes 1979; Rybicki & 
Hummer 1991; Choi et al. 1995; Juvela 1997; Park & Hong 
1998; Hogerheijde & van der Tak 2000). 

In this paper observations are presented of the proto- 
stellar candidates NGC1333-IRAS2 and Serpens SMM4, in 
transitions of HCO+, H1^C0+, CS, CO, "CO and C^O. 
The HCO^ and CS transitions preferentially trace high den- 
sity gas, whereas CO traces a much wider range of gas den- 
sities. The paper is laid out as follows: Section 2 introduces 
the A-iteration method of numerical radiative transfer; Sec- 
tion 3 describes our approach to the modelling; Section 4 
explores the sensitivity of the model to the various free pa- 
rameters; Section 5 describes our observations and data re- 
duction; Sections 6 & 7 present the results of our observa- 
tions for NGC1333-IRAS2 and Serpens SMM4 respectively; 
Section 8 compares the observations with the model predic- 
tions and finds the best fits to the data; Section 9 presents a 
brief summary of our main findings. The non-expert reader 
may prefer to read the second half of the paper first, starting 
from Section 5. 



2 NUMERICAL RADIATIVE TRANSFER 

Consider a cloud of gas with a specified distribution of den- 
sity, kinetic temperature and composition, which may have 
internal turbulent and systematic motions. Let any radi- 
ation sources not forming part of the cloud itself also be 
specified. For each molecular species, there exists a steady 
state solution for the distribution of rotational energy level 
populations as a function of position in the cloud. In the 
modelling of this paper we numerically calculate this dis- 
tribution for the idealised case of a spherically symmetric 
model cloud, and predict the observed line profiles. This is 
a complex problem, due to the fact that well separated parts 
of the cloud can interact radiatively with each other. The A- 
iteration method, described below, is conceptually one of the 
simplest techniques for solving this kind of problem. 



2.1 The A-iteration method 

The method is begun by choosing an initial radiation field 
in a more or less arbitrary manner. From this a 'false' run 
of level populations may be calculated, by assuming the va- 
lidity of the steady state rate equations. Radiative transi- 
tions between these level populations will generally produce 
a radiation field which departs from the one originally as- 
sumed. If this radiation field is used to calculate a new set 
of level populations, and the procedure is repeated a suffi- 
ciently large number of times, the radiation field and level 
populations should eventually converge on a mutually con- 
sistent solution. A-iteration is a kind of diffusion process, 
where imbalances in the radiation field are smoothed out 
over a length-scale corresponding to approximately one op- 
tical depth at each iteration step. 



A number of radiative transfer models have been pub- 
lished (e.g. Rawlings et al. 1992; Zhou 1992; Walker et al. 
1994). However, the modelling in this paper was carried out 
using a modified version of the Stenholm A-iteration code 
developed by Stenholm and subsequently expanded by Little 
and co-workers at the University of Kent (Stenholm 1977; 
Matthews 1986; Heaton et al. 1993). The code uses the above 
method to solve the spectral line radiative transfer problem 
for the rotational transitions of linear molecules in a spher- 
ically symmetric model cloud. Radial profiles of systematic 
velocity, temperature, density, tracer molecule abundance 
and micro-turbulent velocity dispersion may be specified. 

The model cloud is discretised using a number of spher- 
ical shells, and the level populations in each shell are de- 
termined in the iteration process from the calculated mean 
radiation intensity in the co-moving shell frame. Once the 
level populations have converged, a calculation is made of a 
simulated spectral line observation on the cloud. This part 
of the code carries out straightforward numerical integra- 
tions of the equation of transfer along parallel lines of sight 
through the cloud to find the emergent intensities for a grid 
of impact parameters, which are then weighted according to 
a two-dimensional gaussian beam profile function. 

2.2 Testing for convergence 

It is important to find a reliable test for the convergence 
of the model. One possible condition is that the maximum 
fractional change rj{i, J)jv, of the level population n{i, J) be- 
tween iterations A'' — 1 and A'', is less than some specified 
value q: 



,. .,, n(i, J)]v - n(i, J)jv_i 



(1) 



n(i, J)jv 

for all i and J, where the subscripts specify the iteration 
numbers. In many cases, the value of rj{i, J) is a good esti- 
mate of the absolute fractional error in the corresponding 
level population, in which case convergence is ensured by 
setting a suitably low value of a (a value of 0.02 or less is 
usually found to be sufficient). However, when the rate of 
convergence is particularly slow, this test may give a mis- 
leading result, since the fractional change between two suc- 
cessive iterations may then significantly underestimate the 
total change produced over a large number of subsequent it- 
erations. The rate of decrease of 77 should somehow be taken 
into account. 

Dickel & Auer (1994) use the following procedure to es- 
timate the total fractional error r;tot(4, J)]v of the level pop- 
ulation n{i, J) after the A'th iteration: 



r/tot(i, J)]v 



r;(j,J)i> 



[»?(«, J)]v/r;(J,J)]v-i] - 1 



(2) 



The convergence criterion is then ?7tot(i,J) < a for all i 
and J, for some specified a much less than 1. Note that if 
ri{i,3)N ^ ri{i,3)N-i then 77tot(«,J)jv is much larger than 
7)(i,J)jv - i.e. the total error is much larger than the frac- 
tional change between successive iterations. 

The rate of convergence is fastest when the optical 
depth of the cloud is small, or when the hydrogen number 
density is much greater than the critical thermalisation den- 
sity. When the optical depth is small, radiation can travel 
across the whole cloud in one iteration, allowing the level 
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Figure 1. Comparison of analytical and numerical predictions of the CS(J= 5 — > 4) spectra from four uniform LTE models with 
progressively increasing CS abundance. The solid and dashed lines show the numerical and analytical predictions respecitively. Each 
model has a radius of 10,000au, a uniform kinetic temperature of 20K, a uniform hydrogen number density of lO^^cm"'', and a uniform 
FWHM velocity width of 1.0km s~^. The uniform CS number density in the first model is 10~^cm~'^, and increases by factors of 10 in 
successive models, to 10~^cm~^ in the last model. The microwave background intensity has been subtracted from both sets of spectra. 



populations in different parts of tfie cloud to adjust to eacli 
other quickly. When the density is close to, or above, the 
critical density for the principal thermally accessible transi- 
tions of the molecule, the level populations are determined 
mainly by molecular collisions, and are insensitive to changes 
in the radiation field. 

For optically thick regions where the density is below 
the critical density, the convergence rate is generally much 
slower. The low density results in low collision rates, and 
the level populations therefore depend more sensitively on 
the radiation field, while the high optical depth means that 
the radiation field will take many iterations to adjust to the 
conditions in distant parts of the cloud. The code converges 
most rapidly near to LTE and when optically thin. 

We have made several alterations to the code since it 
was made available to us. These have included changes to 
the internal structure of the code by, for example, improving 
the shell subdivision algorithm and incorporating rotation 
into the line profile calculation. However, the former did not 
affect the conclusions, it merely gave the code more reso- 
lution. Some of the more significant changes made are now 
discussed in detail. 



2.3 The Rybicki approximation 

An important modification made to the radiative transfer 
calculation was the removal of the approximation of core sat- 
uration used in the original code, known as the Rybicki ap- 
proximation (Rybicki 1972; Stenholm 1977; Rybicki 1984), 
which may be summarised as follows: suppose that at some 
point in the cloud the following relation for the optical depth 



r^ at frequency i' is satisfied for all directions ro and solid 
angles k: 



T^{ro,Zb,u,'k) > 7, 



(3) 



where zt is the distance to the boundary of the model and 7 
is some prescribed constant. The intensity I in all directions 
at that point is defined in the Rybicki approximation to be 
equal to the local source function S, i.e: 



I^{ro,iy) = Si,{ro,u). 



(4) 



This equation is only applied for the particular values 
of f for which t^ > j. For increasing values of 7, the vol- 
ume of the cloud to which the approximation applies be- 
comes smaller, and disappears altogether when 7 increases 
beyond the maximum optical depth of the cloud. For the 
problems we consider, the source function cannot be as- 
sumed to be uniform over distances corresponding to a few 
optical depths, and the Rybicki approximation cannot be 
used. We therefore do not use the Rybicki approximation in 
the modelling in this paper. 



2.4 Solid body rotation 

We have incorporated rotation into the code, in the form of 
'solid body' rotation - i.e. the angular velocity of all points in 
the cloud about some fixed axis is a constant, but systematic 
radial motions in the rotating frame are still allowed. The 
radiative transfer effects of solid body rotation were fully 
accounted for by changes to the section of the code which 
calculates the emergent line profile from the cloud, after the 
A-iteration solution has been found. No change needs to be 
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Figure 2. Comparison of the CS line profile calculations from figure 13 of Kriigel & Chini (1994) — top - with those produced by 
Stenholm for the same model cloud - lower. For each transition, line profiles are plotted for a number of different FWHM beam sizes, 
as indicated in the bottom right hand panel of the top diagram. The model cloud has an outer radius of 6684 au and is assumed to lie 
at a distance of 450 pc. The adopted hydrogen number density and velocity profiles are of the form ngjC^) ocr^^" and v(r) ocr^^" 
respectively, with njjj = 1-1 X 10^ and v=— 0.237kms~^ at the outer edge of the cloud. The kinetic temperature is set to lOK up to 
the half-point radius, rising linearly to a value of 20K at the outer edge of the cloud. The microturbulent velocity width and the CS 
abundance are assumed to be constant throughout the cloud, with values of 0.5kms~^ (FWHM) and 2xl0~® respectively. 



made to the A-iteration code itself, since the solution for the 
level populations in the cloud does not change when solid 
body rotation is introduced, as is now explained. 

The radiation field in a frame co-moving with the gas at 
a given point in the model can only be affected by a super- 
imposed velocity field if the latter introduces changes in the 



relative line of sight velocities of other parts of the cloud (ig- 
noring relativistic effects of order v^ ji? or less). Solid body 
rotation does not introduce changes in relative line of sight 
velocities from one part of the cloud to another. For ex- 
ample, consider a non-rotating cloud which has systematic 
radial velocities Vr(ri) and Vr(r2) at the positions ri and 



© 0000 RAS, MNRAS 000, 000-000 



Modelling protostellar infall 5 



r2 respectively. Writing ri2 =j r2 — ri | for the distance sep- 
arating the two points, the component of relative velocity, 
vios, along the line of sight between ri and r2 is given by: 



vios = (r2-ri) • [vr(r2)-Vr(ri)]/ri2. 



(5) 



Now assume the cloud is made to rotate in solid-body 
fashion with angular velocity ft. Let the new systematic 
velocities at positions ri and r2 be v(ri) and v(r2), where: 



v(ri) = Vr(ri) + Jl X ri, 
v(r2) = Vr(r2) + $7 X r2. 



(6) 

(7) 



Substituting into equation 5, the terms in fl disappear, 
and the final expression is unchanged. This implies that the 
solid body rotation has no effect on the relative line of sight 
velocities between any two points in the cloud, and hence 
the radiation field in the co-moving frame at each point in 
the cloud is not affected by solid body rotation. This result 
only applies to the radiation emitted by the cloud itself. 
Radiation entering the model from outside will be Doppler 
shifted in the gas frame by the rotation, but the effect of 
this is completely negligible for realistic external radiation 
fields. 

However, solid body rotation does make a difference to 
the prediction for the observed line profile, since the rota- 
tion does introduce changes in the line of sight velocities of 
different parts of the cloud relative to an observer outside 
the cloud. If z is the unit vector parallel to the line of sight 
from the observer to the cloud, and Vz (r) is the line of sight 
velocity (relative to the observer) of the gas at the position 
vector r (relative to the centre of the cloud), then: 



Vz(r) = z • Vr(r) + z ■ (fi X r) = Vr,z(r) + yfl^ 






(8) 



where Vr,z(r) is the component along the line of sight of the 
radial component of the gas velocity at r, and x and y are, 
as before, spatial co-ordinates in the plane perpendicular to 
the line of sight at the distance of the cloud, fix and Qy are 
the components of the angular velocity vector in this plane. 
The component of angular velocity parallel to the line of 
sight, fiz, does not appear in equation 8 and has no effect 
on the observed spectrum from the cloud. 



2.5 LTE comparison 

Analytical solutions of the equation of transfer are in general 
unobtainable for all but the simplest of problems. However, 
in the limit that the level populations are completely ther- 
malised by collisions - i.e. in local thermodynamic equilib- 
rium (LTE) - the spectral line calculations from the numer- 
ical model may be compared directly with analytical predic- 
tions. 

In LTE, the excitation temperature for each transition 
approaches the kinetic temperature Tkin of the gas, and the 
source function for the transition J-^ J— 1 becomes S'^(r) = 
By[i^j,rkin(i')l- If Tkin is constaut over the cloud, then the 
solution to the equation of transfer is: 

/4vz,J)=jr(J)e-"-<"-'+B.(j^.T,rkin)(l-e-"-(^^>), (9) 

where t^{vz) is the optical depth of the cloud at the fre- 
quency corresponding to the velocity bin v^- For simplicity, 
we consider the intensity emerging from a uniform density 
isothermal cloud layer of thickness L, with no systematic 



motions, and with a constant velocity dispersion (uv) and to- 
tal tracer molecule number density (ntot). The optical depth 
through the cloud layer may then be written as: 



r,y{vz) = k{u)L 

hujL 

An 
hvjL 



-->r^ 



=>T^ = 



47r 



(10) 
[B,„(J)n(J - 1) - B„KJ)w(J)l'^(v.) (11) 



B„,n(J)[exp(/ii/j/A:rkin) - l]<^(v,), (12) 



where n(J) and n(J — 1) are the number densities of tracer 
molecules in the upper and lower level of the transition, and 
other symbols take their usual meanings. In the last equation 
we have used the fact that in LTE: 



n(J - 1) _ ff(J-l) 

«(J) " ff(J) 
and also: 



exp(/ii/j/fcrkin), 



fl(J)B„,(J)=p(J-l)B,„(J). 

The explicit equation for <J}{yz) is given by: 



^^!/,j\/27rcxp ( ~j^ 



(13) 



(14) 



(15) 



and the number density n( J) of tracer molecules in the level 
J may be evaluated by: 



n(J) = ntot 



(2J + l)exp(-/it/.j/fcrkin) 
E"=o(2J' + 1) exp(-ft^j,/A:Tki„ 



(16) 



The sum in the denominator is the partition function, and 
usually converges rapidly after about ten or so terms (de- 
pending on Tkin), and can be easily evaluated to the required 
accuracy. 

Substituting these expressions into equation 9, with 
/""'(J) = Bv{i^j,2.73K), gives an expression for the pre- 
dicted intensity as a function of Vz- To compare these pre- 
dictions with the Stenholm code, we used a static model 
cloud with a uniform kinetic temperature, hydrogen number 
density, velocity dispersion, and abundance. The hydrogen 
number density was set several orders of magnitude above 
the critical density to enforce LTE. 

Figure 1 shows comparisons of the analytical and nu- 
merical predictions for the CS(J= 5^4) line, from LTE 
models with progressively increasing CS abundance, cov- 
ering a wide range of cloud optical depths. The excellent 
agreement between the predicted spectra in all parts of the 
line, and across the large range of optical depths, is a good 
indication that the integration of the equation of transfer 
along the line of sight is being carried out correctly. 



2.6 Non-LTE comparison 

The LTE comparison shown above, while instructive, does 
not test the ability of the code to carry out non-LTE radia- 
tive transfer calculations reliably. We have not found any ex- 
act analytical solutions for non-LTE radiative transfer prob- 
lems in spherical geometry, to compare with the code in a 
similar manner as shown above. Instead, we attempt to re- 
produce the line profiles obtained in a previous non-LTE ra- 
diative transfer study (Kriigel & Chini 1994, hereafter KC), 
which used a somewhat different numerical method. 

Figure 2 shows a comparison of profiles taken from KC 
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with the corresponding profiles calculated by our Stenholm 
code. The figure caption gives details of the model cloud 
used, which has strong velocity and density gradients to- 
wards the centre, and is far from LTE. Both calculations use 
the CS collision rates of Green & Chapman (1978). Overall, 
the agreement between the two sets of line profiles is very 
good. The intensities of the line profiles predicted by Sten- 
holm are up to 10% weaker than the KC profiles, and the 
disagreement appears to be greatest for the smallest beam 
sizes. 

We checked our calculations, and found that this was 
not the result of inadequate sampling of the lines of sight 
used in the beam convolution, by re- calculating the profiles 
for a range of different beam sampling densities. The small 
disagreement in the line profile strengths is not serious, and 
may be caused by differences in the spatial, angular or fre- 
quency discretisation schemes used in the models, or small 
differences in the precise values of some of the parameters 
used. Particularly encouraging is the very good agreement 
in the shapes of the corresponding line profiles. 



3 APPROACH TO THE MODELLING 

In this section we discuss our approach to the modelling 
of infall line profiles, and investigate the qualitative depen- 
dence of the predicted line profiles on a number of different 
model parameters. The approach to the modelling is neces- 
sarily a compromise between the desire to represent realis- 
tically the object being modelled, and the practical need to 
limit the parameter space and running time of the numer- 
ical simulations. The observations are modelled using the 
spherically symmetric micro-turbulent A-iteration radiative 
transfer code described above. Even with the considerable 
simplification of spherical symmetry the number of model 
parameters is formidable. 

Radial profiles of density, kinetic temperature, molecule 
abundance, systematic velocity, and micro-turbulent veloc- 
ity width must be specified. Additional parameters include 
the choice of inner and outer boundaries of the model, the 
incident radiation fields at the inner and outer boundaries, 
and the distance to the object. Some assumptions must be 
made about the forms taken by these profiles, so that they 
may be described by specifying only a few parameters. Some 
of the profiles are assumed to have a power law dependence 
with radius, which requires only the power law exponent 
and a normalisation factor to be specified. Wherever pos- 
sible, the choice of parametric forms used to describe the 
radial profiles is guided by theoretical considerations. We 
now discuss the profile of each physical parameter in turn. 



The infall velocity profiles in non-rotating clouds which col- 
lapse from different initial conditions should approximately 
converge, even when the mass accretion rates are very dif- 
ferent, as long as the comparison is made at times when the 
central protostellar masses are equal (Hunter 1977). Hence 
we use the singular isothermal sphere (SIS) model velocity 
profiles to parameterise the infall velocity field in the radia- 
tive transfer modelling, which allows the velocity field to be 
described in terms of only two quantities: the effective sound 
speed ttcff and the infall radius rinf. We stress that these are 
merely used as parameters characterising the velocity field, 
and their actual physical significance is model dependent. 

We have found the following empirical fit to the exact 
velocity profile predicted by the standard SIS model, for 



Vr{r) 

flcff 



-^/2 



r 

Hnf 



+ \/2 



r 

Hnf 



(18) 



This equation fits the SIS collapse field to within a few per 
cent. If the collapsing cloud has some initial rotation, then 
centrifugal support may become important at small radii, 
leading to a decrease in the radial velocity close to the centre 
of collapse. With an initial angular velocity O, then the inner 
radius, ra, at which centrifugal forces begin to dominate 
depends on the mass already accreted, Mace, and the radius 
R enclosing this mass (Terebey, Shu & Cassen 1984) through 
the following relation: 



Td = 






(19) 



3.2 Micro-turbulent velocity profile 

As discussed above, the radiative transfer code used to 
model the data in this paper makes use of the micro- 
turbulent approximation. This assumes that random 'tur- 
bulent' motions can be treated in the same way as ther- 
mal molecular motions, by incorporating them into the local 
line profile function. Comparisons between radiative transfer 
codes which model macroscopic and microscopic turbulence 
suggest that the main effect of relaxing the micro-turbulent 
approximation is to reduce the strength of self-absorption 
features (e.g. Park & Hong 1995). 

Although observations appear to indicate random mo- 
tions of some kind must be present in protostellar envelopes, 
very little theoretical work has been done on the role of 
turbulence and turbulent support in protostellar collapse. 
Lizano & Shu (1989) derived a 'turbulent equation of state' 
using the well known 'Larson's Laws' (Larson 1981), which 
relate the observed velocity dispersion (a), size (R), and av- 
erage density (p) in molecular clouds and cores: 



3.1 The systematic velocity profile 

Following the onset of fully dynamical collapse, the radial 
velocity profile should tend towards a freefall profile (e.g. 
Larson 1969; Shu 1977; Hunter 1977; Zhou 1992; Foster & 
Chevalier 1993). The pure free-fall velocity profile may be 
written as: 



Vr{r) 



{^y 



(17) 



R-\ 



(20) 
(21) 



These correlations apply to ensembles of clouds and also 
to observations of individual clouds in different molecular 
tracers, and are associated with the turbulent motions in 
the clouds (e.g. Myers 1983). Combining the two relations 
above gives: 



a (X p 



(22) 
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Zhou (1992) used this relation to infer the variation of tur- 
bulent velocity with density in his radiative transfer calcu- 
lations of infall line profiles. 

Vazquez-Semadeni, Canto & Lizano (1998) have re- 
cently argued that Larson's Laws may be inappropriate 
to dynamically collapsing clouds, because they were estab- 
lished from observations of clouds which were in equilibrium 
configurations. These authors carried out numerical hydro- 
dynamic and magneto-hydrodynamic simulations to test the 
behaviour of the turbulence during the collapse, and found 
that, in contrast to equation 22, the turbulent velocity dis- 
persion actually increases with density during the collapse. 

For purely hydrodynamic and weakly magnetic collapse 
models the relation between the turbulent velocity disper- 
sion and the density was seen to approach a power law of 
the form a oc p^'^ (i.e. Ptb oc p^). In the case of strongly 
magnetically inhibited collapse, the power law index was ob- 
served to change to (t ex p^'^ (Ptb oc p^ ), consistent with 
the predicted behaviour for the slow compression of Alfven 
waves (McKee & Zweibel 1995). The behaviour of the turbu- 
lence in infalling protostellar envelopes is thus theoretically 
rather uncertain. Turbulence nevertheless plays a very im- 
portant part in the formation of line profiles, and is one of 
the principal sources of uncertainty in the infall line profile 
modelling. 



3.3 Density profile 

In a region where the infall velocity has a free-fall profile 
there is often a strong tendency for the density profile to 
adopt an r^^" distribution, since if the rate of injection of 
material into the accretion stream is not changing rapidly 
with time, the quantity Aivr^pVr should be approximately 
independent of radius. The prediction for the density profile 
deep inside the infalling region may be written as: 



p{r) = 



7rV32 Gr 



1/2 



(23) 



where, as before, rinf = acst is the infall radius, and Oeff is 
the effective sound speed. The exact model density profile 
inside the whole infalling region is very well approximated 
(within ~1 per cent) by the following empirical formula: 



p{r) = 



2nGrl, 



0.35 



rinf 



0.65 



r-inf 



(24) 



Outside the infall radius, the density profile is given by: 



p{r) 



27rGr2 ' 



(25) 



There are several competing theoretical models which pre- 
dict different forms of the density profiles at various stages 
in the contraction and collapse of a dense cloud core (e.g. 
Foster & Chevalier 1993; Basu & Mouschovias 1994; Whit- 
worth & Ward-Thompson 2001). 

However, once a central protostar has formed there is 
a tendency in the theoretical models for the forms of the 
density and velocity profiles to converge towards the SIS 
model forms, in particular towards the centre of the collapse 
(e.g. Foster & Chevalier 1993). We therefore initially search 
for fits to our observations using density profiles of the above 
form, although we do not require that the adopted density 



profile uses the same values for the parameters Tint and Oeff 
as used in the velocity profile. 



3.4 Kinetic temperature profile 

The temperature of the gas in an infalling protostellar en- 
velope is determined by the balance of heating and cooling 
processes. Ceccarelli, HoUenbach & Tielens (1996) modelled 
these processes to derive the temperature profile in pro- 
tostellar envelopes with central embedded heating sources. 
They found that in all the cases studied, the gas temperature 
profile stays within ~30% of the dust temperature, usually 
lying slightly below the dust temperature, Tj. In the follow- 
ing it is assumed that the gas temperature is well coupled 
to the dust temperature, although this assumption may be 
significantly in error for hydrogen number densities below 
~ lO'^cm-^ 

If the emitted radiation reaching the optically thin part 
of the envelope peaks at a wavelength A*, then the temper- 
ature of the dust at a radius r will be given by: 



Td(r) ex 






l/(4+/3) 



(26) 



The value of A* depends on the details of the dust radia- 
tive transfer through the optically thick inner region (i.e. 
how the material in this region is distributed). Kenyon et 
al. (1993) used the diffusion approximation in the optically 
thick region to estimate the value of A,t. For a density profile 
of the form p{r) — por~' this gives: 



A* ex 



Po 



2/(9-1) 



1/(4+4/3) 



(27) 



Substituting this expression into the above equation and set- 
ting P — 1.5 and q = 1.5 (appropriate for a hydrodynami- 
cally collapsing region), we find: 



Td{r) oc po 



0.055 r0.21^ 



(28) 



The temperature profile in the optically thin part of the 
envelope is therefore only weakly dependent on the density 
normalisation and the luminosity (see Kenyon et al. 1993 for 
further discussion). 

For the objects modelled in this paper, the optically 
thick region is expected to be very much smaller than the 
beam size of the observations and is comparable to the inner 
shell radius used in the radiative transfer models. At suffi- 
ciently large radii, the gas temperature will approach the 
ambient temperature of the larger scale cloud. We therefore 
adopt the temperature profile appropriate to the optically 
thin limit (T oc r" where —0.4 < s < —0.33) in the inner 
region, and switch to a flat profile at the radius where the 
power law inner profile falls below the ambient cloud tem- 
perature. 



3.5 Abundance profile 

The profiles of CS and HCO^ relative abundance in an in- 
falling protostellar envelope are governed by time depen- 
dent chemistry and the transfer of molecules between the 
gas phase and the ice mantles surrounding dust grains. This 
problem has been studied by Rawlings et al. (1992) and 
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Figure 3. Predicted JCMT line profiles for the canonical infall model (see text for details) for an assumed distance of 200 pc. The rare 
isotoponier lines are plotted as dashed lines in the top panels. 



HCO*{J = 4-3) 



;.. 



/v 



f.'-\ 

! --i i \ 

! I \ I \ 







V(km/s) 



CS(J = 5-4) 




V{km/s) 



Figure 4. Plots showing the contribution of the emission from the envelope hemisphere nearer to the observer (dashed line) and the 
attenuated emission from the further half of the envelope (dotted line), to the predicted HCO"'"(J = 4^3) and CS(J = 5 — > 4) line 
profiles (solid lines). 
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Figure 5. Plots showing the unattenuated emission from the further half (from the observer) of the model envelope (dot-dashed line) 
in the HCO"*'(J = 4 — > 3) and CS(J = 5^4) lines. The full line profile (solid line) and the attenuated emission from the further half 
(dotted line) are shown for comparison. 
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Willacy, Rawlings & Williams (1994) in the case of isother- 
mal collapsing clouds. The objects studied in this paper con- 
tain internal stellar heating sources, and are unlikely to be 
well described by an isothermal model. The chemistry of 
internally heated infalling protostellar envelopes has been 
studied by Ceccarelli, HoUenbach & Tielens (1996). We re- 
fer the reader to all of these references for further discussion 
of the way in which chemistry can affect observed abundance 
profiles. 

Given these uncertainties and the requirement to limit 
the parameter space, we used flat abundance profiles in our 
modelling. The fractional abundances of CS and HCO"'" 
measured in dense molecular gas cores occupy similar 
ranges, typically between 10~^ and 10^* (e.g. Irvine et al 
1987; McMullin, Mundy & Blake 1994; Blake et al. 1995). 
The measured values of [^^C/^^C] at the galacto-centric ra- 
dius of the sun (which applies to all of our observed objects) 
lie in the range 50-100 (Wilson & Rood 1994). We there- 
fore restrict the abundances used in the modelling to these 
ranges of values. 



4 PARAMETER STUDY OF INFALL 
PROFILES 

4.1 The canonical infall model 

It is useful to define a specific infalling envelope model, with 
plausible physical parameters, which may be used to inves- 
tigate some general properties of infall spectral line pro- 
files. We use the above density and velocity profiles at a 
stage when a mass M0/2 has already accreted onto the cen- 
tral protostar, and a further ~M0/2 of envelope gas is in- 
falling towards it. If we choose an effective sound speed of 
fflcff ~ 0.35 kms~^, then an infall radius of ~3700 AU is im- 
plied. We therefore use the following model relations for the 
systemic velocity and hydrogen number density profiles: 



Vr = 0.49 



n/fo 



(3705I!7)""-°-^K 



= 0.35 



3700AU 



3700^17 



+ 0.65 



km/s. 



S700AU 



2 X 10 Vcm3 

and outside the infall radius the systematic velocity is set to 
zero. The density profile is given by: 



nH2(r) = 2.0 X lO' 



V3700AUJ 



(29) 



.3700AUy 

We truncate the density profile at an outer radius of 
10000 AU, which encloses a total mass of 2.75 M©. This cor- 
responds to the typical radius at which dense cores merge 
with the more diffuse gas in the ambient molecular cloud. 

The adopted micro-turbulent velocity dispersion is 
o"tb = 0.3 km s^^, which is assumed to be spatially constant. 
Although there is little empirical or physical justification for 
a flat turbulent velocity profile (see the above discussion), 
we make this assumption as a zeroth order approximation, 
in the absence of any clear consensus as to the correct de- 
scription of turbulence in infalling envelopes. 

As discussed above, the kinetic temperature profile (as- 
sumed to be equal to the dust temperature profile) in the 
optically thin part of the envelope is expected to be fairly 
insensitive to the luminosity of the central source. We choose 
a canonical temperature profile of: 



Molecule 


Mass 
(10-27kg) 


(10-30Cm) 


Brot 

(GHz) 


hBrot/k 
(K) 


HCO+ 

H13CO+ 

CS 


48.16 
49.81 
73.02 


12.2 
13.0 
6.54 


44.59 
43.38 
24.50 


2.14 
2.08 
3.50 



Table 1. Assumed physical constants of the molecules modelled. 



T{r) = 20 ( 



lOOOAU 



K. 



(30) 



The normalisation of the temperature profile was fixed us- 
ing the dust temperature profile calculated by Kenyon et 
al. (1993), for a source luminosity of ~ 151/q. The varying 
kinetic temperature in the envelope is not strictly consis- 
tent with the assumption of constant aog implied by the SIS 
model, where a^g — a^^^ + kT/fh. However, our primary con- 
cern is not to investigate the SIS model solution in detail, 
but rather to investigate the essential features of line forma- 
tion in infalling envelopes, whether the infall is described by 
the SIS model or not. 

We set both the CS and HCO^ relative abundances to 
5.0 X 10"^ and assume a [H^2CO+]/[H^^CO+] isotopic ratio 
of 60. The set of molecular constants used in the radiative 
transfer modelling are listed in Table 1, where /i denotes 
the permanent molecular dipole moment and other symbols 
take their usual meanings. The data are taken from the JPL 
spectral line catalogue (Poynter & Pickett 1985). 



4.2 The origin of the infall profile 

We now use the radiative transfer code described above to 
simulate observed spectral line observations of the canonical 
infall model, and investigate how the predicted line profiles 
depend on a number of different model parameters. Figure 3 
shows the predicted HCO+ , H^^CO+ and CS line profiles 
for the canonical infall model described above. The model 
cloud was assumed to lie at a distance of 200 pc, and the 
appropriate beam sizes for each transition were used when 
carrying out the beam convolution. 

Double-peaked, blue-asymmetric line profiles are seen 
in all of the main isotopomer transitions, apart from the 
CS(J = 7-^6) line, which is nevertheless skewed slightly 
bluewards of the systemic velocity. The minima between the 
peaks of the double-peaked profiles lie close to the systemic 
velocity. The HCO^(J = 3^2) line shows a stronger self- 
reversal than the HCO^(J = 4^3) line, and in both of 
these lines the self-reversal is stronger than the CS(J = 5 ^ 
4) line. The two single-peaked H^'^CO^ lines peak very close 
to the systemic velocity, coinciding with the minima in the 
main line profiles. Low level high-velocity wings are visible 
in all of the main line profiles, and appear to be stronger in 
the higher rotational transitions of both molecules. 

To illustrate how the asymmetries in line profiles arise, 
we show in Figure 4 the separate contributions of the near 
and far portions of the envelope to the line profiles for the 
HCO+(J = 4^3) and CS(J = 5^4) transitions. 
The nearer hemisphere (to the observer) of the envelope (in 
which the systematic velocity is red-shifted) is seen to be 
mostly responsible for the red-shifted peak in the line pro- 
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Figure 6. Plots showing the contribution to the total beam-averaged optical depth through the envelope (solid line) from the near and 
far hemispheres of the canonical infall model (dashed and dotted lines respectively), in the HCO+(J = 4^3) and CS(J = 5 — » 4) lines. 
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Figure 7. Plots illustrating how the HCO+(J = 4 — > 3) evolves as progressively more of the cloud is included in the line profile 
calculation. The profiles represent different 'slices' through the model cloud, where the slice planes are perpendicular to the line of sight 
and only the part of the cloud which lies beyond the slice plane is included. The left hand plot shows profiles obtained for slice planes at 
various positions in the far half of the cloud. The distances of the slice planes from the cloud centre along the line of sight are: 5000 AU 
(solid line); 2000 AU (dashed line); 1000 AU (dotted line); and OAU (dot-dashed line). The right hand plot shows the profiles obtained 
for slice planes lying in the near half of the cloud. The distances of the planes from the cloud centre are OAU (solid line); 1000 AU 
(dot-dashed line); 2000 AU (dotted line); 5000 AU (dashed line); and 10000 AU (double-peaked solid line). The radius of the cloud is 
10000 AU. 



file, as expected. However, the more surprising observation 
is that the nearer hemisphere also dominates the emission 
at velocities well into the blue-shifted side of the line profile, 
and contributes significantly to the blue-shifted peak. 

Emission from the far hemisphere only becomes domi- 
nant in the most blue-shifted part of the line. It is interesting 
to note that the contribution of the far hemisphere to the 
blue peak is no greater than the contribution of the near 
hemisphere to the red peak, and it is the extra contribution 
of the near hemisphere to the blue-shifted emission which 
produces the blue asymmetry. This contrasts with the ex- 
planation usually given for the infall asymmetry, in which 
it is argued that the optical depth to the warm dense blue- 
shifted gas just beyond the central protostar is less than the 
optical depth to the red-shifted gas just in front of the pro- 
tostar, due to the Doppler shift of the intervening infalling 
envelope. If this explanation were the correct one, then we 
would expect the blue-shifted emission from the far hemi- 
sphere to be, by itself, much brighter than the red-shifted 
emission from the near hemisphere, which contradicts our 
findings. 

Figure 5 shows how the emission from the far half of the 
envelope is modified by the optical depth of the near half 
of the envelope on its way to the observer. Self-absorption 
clearly plays a very important role in shaping the emer- 
gent line profile. In Figure 6 the contribution of the beam- 
averaged optical depths through both hemispheres of the en- 



velope for both CS(J = 5^4) and HCO+(J = 4^3) are 
plotted. The total optical depth of the HCO+(J = 4^3) 
line is much greater than the CS(J = 5^4) line, which 
accounts for the stronger self-absorption features in the 
HCO^ line. As expected, there is a separation in velocity of 
the peaks of the optical depth profiles from the near and far 
hemispheres, although this is rather small (~ 0.2 kms^^). 

To shed further light on how the asymmetric infall pro- 
files are formed, we show in Figure 7 how the HCO^(J — 
4^3) line profile changes as successive 'layers' of the model 
envelope are included, from the far edge of the envelope for- 
wards. For all the 'slice planes' in the rear of the envelope 
(left hand panel of Figure 7) the predicted profiles are single 
peaked. This is explained by the fact that in the far half of 
the envelope, the higher excitation gas always lies in front of 
the lower excitation temperature gas. Figure 8 shows the ra- 
dial dependence of the systematic velocity, kinetic tempera- 
ture, and HCO+(J = 4 -> 3) and CS(J = 5^4) excitation 
temperature in the model. The fact that the excitation tem- 
perature lies well below the kinetic temperature over most 
of the envelope is a demonstration of the inapplicability of 
the LTE approximation to this problem. 

As successive layers of the front half of the envelope are 
included in the line profile calculation (right hand panel of 
Figure 7), there is a reduction in the intensity in the core 
of the line. This is due to the fact that in the near half of 
the cloud, lower excitation temperature gas lies in the fore- 
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Figure 8. The left hand panel shows the radial profile of systeinatie velocity for the eanonical infall model. In the right hand panel, 
the solution for the exeitation temperatures of the HCO+{J = 4 — > 3) (dashed line) and CS(J = 5^4) (dotted line) transitions are 
plotted. The kinetic temperature profile (solid line) is also plotted for comparison. 
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Figure 9. Plot showing the dependence of the line profiles on the impact parameter of the beam direction relative to the centre of the 
cloud. The solid lines show the on-source line profiles. The line profiles for impact parameters of 1000 AU (dotted line), 2000 AU (dashed 
line), and 4000 AU (dot-dashed line) are also plotted. The assumed distance to the cloud is 200 pc. 
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Figure 10. Plots showing the dependence of the predicted line profiles on the assumed distance to the cloud. Spectra are plotted for 
distances of 50 pc (solid lines), 100 pc (dashed lines), 200 pc (dotted lines), and 400 pc (dot-dashed lines). 
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Figure 11. Plots illustrating the dependence of the line profiles on infall velocity. The infall velocity profiles were calculated for infall 
radii of 2000 AU (solid lines), 4000 AU (dashed lines), 6000 AU (dotted lines) and 8000 AU (dot-dashed lines). The canonical values were 
used for all other model parameters. 
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ground, and absorbs the emission from the higher excitation 
gas behind it. Nevertheless, the emission of this gas is also 
very important in determining the form of the emergent line 
profile. The emission from gas near the centre, in the near 
half of the envelope, is seen in Figure 7 to produce an in- 
crease in the line intensity at strongly red-shifted velocities. 
The contribution of this gas to the emission at lower veloc- 
ities is obscured in the figure by the dominating effect of 
absorption towards the line centre, but Figure 7 shows that 
this emission is significant even on the blue-shifted side of 
the line. In the outermost region of the envelope, where the 
excitation temperature is lowest, the gas is predominantly 
absorbing, and simply deepens the self-reversal of the line 
profile. 



4.4 Infall velocity 

Figure 11 shows line profiles for a number of different infall 
radii, in which only the velocity profile is varied, as dis- 
cussed in section 3.1, and the density profile remains fixed 
as for the canonical model. For increasing values of the in- 
fall radius, larger fractions of the envelope take part in the 
infall motion, and the collapse speed at a given radius in- 
creases approximately as r^^^ . The line profiles become more 
strongly skewed towards blue velocities with increasing in- 
fall velocity, mainly as a result of the diminishing intensity of 
the red-shifted peak. The actual velocities of the two peaks 
show remarkably little variation over the range of infall ve- 
locities considered. Predictably, the models with the highest 
infall velocities produce the strongest emission in the line 
wings. In actual observations this infall signature may often 
be obscured by emission from the outfiow. 



4.3 Impact parameter and beam size 

Figure 9 illustrates how the predicted HCO^(J = 4^3) 
and CS(J = 5^4) line profiles depend on the impact pa- 
rameter of the line of sight. As expected, the strength of 
the lines diminishes with increasing impact parameter, as 
the beam samples increasingly less dense and lower tem- 
perature gas. The degree of asymmetry in the line profiles 
decreases with increasing impact parameter, partly because 
of the lower infall velocities at larger radii, and partly be- 
cause of the lower column densities, and hence optical depth, 
along the line of sight. The lines also become narrower with 
increasing impact parameter, as a result of the decreasing 
infall velocity with radius. 

The distances to many nearby protostellar objects are 
often only known to an accuracy of ~50%, so it is important 
to examine how sensitively the predicted line profiles depend 
on the assumed distance. There is of course a degeneracy 
between the assumed distance to the cloud and the beam 
size. 

In Figure 10 we plot HCO+(J = 4 -> 3) and CS(J = 
5^4) line profiles for a number of different assumed 
distances, whilst keeping the beam sizes fixed. The peak 
line temperatures increase approximately as the inverse of 
the assumed distance (or beam size). The increase in the 
strength of the line wings depends approximately on the in- 
verse square of the assumed distance, which suggests that 
the line wing fiux originates from a region on the source 
much smaller than the beam size. This is the explained by 
the fact that in the velocity profile we have adopted, the 
highest infall velocities lie at the smallest radii, and there- 
fore the high velocity wings of the line profile are formed 
in a very small region around the centre of the cloud. The 
different dependence of the line core and line wings on dis- 
tance produces increasingly broad profiles with decreasing 
distance and/or beam size. 

We conclude that distance uncertainties may contribute 
significantly to the overall uncertainty in the radiative trans- 
fer modelling of protostellar envelopes. Gregersen et al. 
(1997) analysed HCO^ observations of a large sample of 
Class sources using a model cloud at a single 'represen- 
tative distance'. This approach will not give reliable results 
when the scatter in the distances of the objects in the sample 
is large. 



4.5 Tracer molecule abundance 

Figure 12 shows how the predicted line profiles vary with the 
assumed relative abundance of CS and HCO"'' in the enve- 
lope. The relative abundance of the tracer molecule largely 
determines the optical depth through the envelope in the 
observed line, and is therefore a very important parameter 
in deciding the overall appearance of the line profile. This is 
apparent from the figure, where the intensities and shapes of 
the line profiles are seen to vary significantly over the factor 
of 10 range of relative abundances covered. As expected, the 
lowest values of the relative abundance (and hence optical 
depth) produce the weakest infall profiles. The progression 
of line shapes with increasing optical depth is from: a single- 
peaked gaussian; to a single-peaked blue-skewed profile; to 
a profile with a blue-shifted peak and red-shifted 'shoulder' 
(or 'red knee'); to a double-peaked profile with a stronger 
blue-shifted peak. As the abundance and optical depth in- 
crease further, the absorption trough deepens, the line peaks 
become stronger and their separation increases. 



4.6 Turbulence 

Figure 13 shows how the line profiles depend on the mag- 
nitude of the turbulent velocity dispersion, oth, when it is 
assumed to be uniform throughout the envelope. As the tur- 
bulent velocity dispersion increases, the most apparent effect 
on the line profiles is to increase the velocity separation be- 
tween the two peaks, as a result of the broadening of the 
absorption profile of the foreground envelope. The total in- 
tegrated flux in the line tends to increase with the turbulent 
velocity dispersion. This is because the peak optical depth 
of the foreground gas is reduced as the molecules are spread 
over a wider range of velocities, and the similar reduction in 
the optical depth of the strongly emitting gas in the centre 
of the cloud allows more of the emission to 'escape'. 

Comparing the HCO+(J = 4^3) with the CS(J = 
5^4) proflles, we see that the effect of changing the value 
of the turbulent velocity dispersion depends strongly on the 
peak optical depth in the transition. When the optical depth 
is small, increasing the turbulence tends to 'smear out' the 
double-peaked profile, whereas for very optically thick lines 
the double-peaked structure remains just as promiment as 
the velocity gap between the two peaks increases. 
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Figure 12. Plots showing the dependence of the predicted line profile on the relative abundance of the tracer molecule. Spectra are 
plotted for HCO+ and CS abundances (relative to molecular hydrogen) of 10^® (solid lines), 2 X 10^^ (dashed lines), 5 X lO^'* (dotted 
lines), and 10~* (dot-dashed lines). 
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Figure 13. Plots showing the dependence of the predicted line profiles on the magnitude of the turbulent velocity dispersion, assumed 
in this case to be constant throughout the envelope. Spectra are plotted for turbulent velocity dispersions of 0.1 km s^^ (solid lines), 
0.2 km s~^ (dashed lines), 0.3 km s~^ (dotted lines) and 0.4 km s~^ (dot-dashed lines). The equivalent FWHM turbulent velocity widths 
are found by multiplying these numbers by 2.35. 
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Figure 14. Plots showing the predicted line profiles for different forms of the micro-turbulent velocity profile. The adopted relations 
between the turbulent velocity dispersion, crti,, and the gas density, p, that were used are: crti, =constant (solid line); a^^ oc p^^ " (dashed 
line); a^^, oc p^'** (dotted line); and (Jtb oc p^ " (dot-dashed line). Each profile was normalised such that the FWHM turbulent velocity 
width (= 2.35(Tti3) at the half-radius of the cloud (5000 AU) has a value of 0.3 km s^^. 
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Note that this finding is in apparent contradiction with 
some previous models. For example, Myers et al. (1996) used 
a simple 'two-slab' geometry to model infalling protostars, 
and claimed that as a protostar evolves, the fraction of the 
envelope that is infalling increases (under the SIS inside-out 
collapse assumptions) , and that this causes the gap between 
the two peaks to decrease, and the absorption trough to 
become narrower. Our modelling indicates that the amount 
of turbulence is what actually determines the gap between 
the peaks. The two scenarios could only be reconciled if tur- 
bulence were to decrease as a protostar evolves. This may 
seem unlikely, since the associated outflow is more likely to 
increase the turbulence as it evolves. However, a recent claim 
has been made (Jayawardhana, Hartmann and Calvet 2001 
- hereafter JHC) that in fact Class sources may be formed 
preferentially in regions of higher turbulence than Class I 
sources. This is an interesting idea that is consistent with 
our modelling, although the other ideas of JHC that a sig- 
nificant fraction of a Class envelope is not infalling, but 
rather static, is not consistent with our modelling of specific 
Class sources (see below). 

Figure 14 illustrates how the line profiles vary with the 
exponent in the assumed power-law relation between the 
turbulent velocity and the density. Guided by the discussion 
on turbulence in the section 3.2, we investigate turbulent ve- 
locity laws of the form crtb oc p^'* and (Jtb oc p^ " (as found in 
the numerical simulations of Vazquez- Semadeni et al. 1998), 
and (Jtb oc p^^" (derived using Larson's Laws). For com- 
parison we also plot the profile produced by the canonical 
model, which assumes a flat turbulent velocity profile. 

The profiles are clearly quite sensitive to the form of the 
turbulent velocity law. The profiles obtained using the crtb oc 
p^^" law (dashed lines) are appreciably narrower than the 
other profiles. This is particularly evident in the centres of 
the H"CO+(J = 4-^3) and CS(J = 7^6) lines. This 
arises because the p~^" law assigns the lowest turbulent 
velocity dispersion to the hot, dense, inner regions of the 
cloud, where most of the line emission takes place. These 
profiles also tend to show the strongest emission close to the 
line centre, mainly arising from the lower optical depth of 
the outer envelope due to the enhanced velocity dispersion 
there. However, in the case of the H"C0"''(J = 4-^3) 
and possibly the CS(J = 7^6) lines, this is due to the 
increased optical depth in the centre of the cloud, as these 
lines are not completely optically thick. 

As the exponent s in the relation crtb oc p" increases 
from —0.5 to -1-0.5, the line profiles become much more 
dominated by broad line wings. For the positive exponents, 
the largest velocity dispersions are located in the centre of 
the cloud, which causes the emission from this region to 
be distributed over a wide velocity range. The low turbu- 
lent and systematic velocities in the outer envelope con- 
fine the absorption of the gas in this region to line cen- 
tre velocities, allowing the enhanced high velocity emis- 
sion from the core to remain virtually unobscured. Despite 
the change in overall appearance, the HCO^(J = 4 — > 3) 
and CS(J = 5^4) profiles retain the characteristic blue- 
asymmetric, double-peaked structure. The strong wings of 
these lines are similar to those normally associated with out- 
flow emission. The more optically thin H^"^CO^(J = 4^3) 
and CS( J = 7^6) lines may provide some means of distin- 
guishing between these two possibilities. If the wings in the 



main line proflles are caused by enhanced turbulence near 
the centre of the cloud, then the optically thin lines sould 
also show very broad profiles. Conversely, if the wings are 
tracing a bipolar outfiow, then the optically thin line pro- 
files, which are expected to be dominated by the envelope, 
should be narrower. 



4.7 Solid-body rotation 

Figure 15 shows how the predicted line profiles change when 
a solid-body rotation is superposed on the infall velocity 
field. This velocity field may not be physically realistic, since 
angular momentum conservation will tend to cause the an- 
gular velocity of the gas to increase as it falls inwards, even 
in the presence of magnetic braking (e.g. Mouschovias 1994). 
We use this example simply to illustrate some of the qualita- 
tive effects of rotation on line profiles. A proper treatment of 
differential rotation would require a three-dimensional radia- 
tive transfer analysis, which is computationally a very large 
step from the spherically symmetric problem presented here. 

The on-source line profiles are least affected by rotation, 
only showing a small decrease in the height of the peaks, and 
a very slight broadening of the line as a whole. For differ- 
ential Keplerian (fisini ~ r~^") rotation, the on-source 
profile might be more strongly affected, since the largest ro- 
tational velocities would then lie closest to the centre. 

The off-centre line profiles are more affected by solid- 
body rotation than the on-source profiles. As expected, the 
shifts in the centroid velocity of the line profiles follow the 
rotational velocity gradient. As well as shifting the centroid, 
the rotation also significantly distorts the shape of the off- 
centre line profiles. On the blue-shifted side of the rota- 
tional velocity gradient, the rotation produces an enhanced 
blue asymmetry in the line profiles, in addition to the blue- 
asymmetry produced by the infall. At the positions on the 
red-shifted side of the rotational velocity gradient, the blue 
infall asymmetry is completely reversed by the rotation. 

Figures 16 & 17 show HCO+(J = 4^3) and CS(J = 
5-^4) centroid velocity contour plots from the canonical 
infall model for projected rotational angular velocities of 
ilsin j=0 and 30kms~^ pc~^ respectively. The non-rotating, 
infalling model produces circularly symmetric negative cen- 
troid velocity contours, which reach a minimum at the cen- 
tral position. The solid body rotation tends to produce ap- 
proximately parallel contours of centroid velocity, aligned 
with the axis of rotation. 

However, the effect of the infall is to cause the con- 
tours to distort towards bluer velocities close to the centre 
of the cloud, producing the appearance a 'blue bulge' in the 
centroid velocity contour plot, which encroaches into the 
red-shifted half of the rotational velocity gradient (Walker, 
Narayanan & Boss 1994; Zhou 1995) 

We can use measurements of the centroid velocity gradi- 
ent across our spectral line maps of protostellar envelopes to 
place limits on the rotation rate of these envelopes. In Fig- 
ure 18 the calculated HCO+(J = 4^3) and CS(J = 5 -> 4) 
centroid velocities are plotted, along axes perpendicular to 
and parallel to the rotation axis, at 1000 AU intervals (i.e. 5 
arcsec intervals for a distance of 200 pc). A straight line fit to 
the centroid velocity along the perpendicular axis is shown 
in the left hand panel of the figure. It is apparent that the 
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Figure 15. Plot showing tlie effect of solid-body rotation on infall line profiles. HCO"'"(J = 4 — > 3) and CS(J = 5^4) line profiles are 
plotted in the left and right hand columns respectively. The solid lines show the predicted line profiles for the canonical infall model with 
a projected angular velocity of solid-body rotation of f2sini=30kms^^ pc^^, where i is the inclination angle of the rotation axis to the 
line of sight. The lines are plotted for several impact parameters. Displacements from the centre of the cloud along a line perpendicular 
to the rotation axis in the plane of the sky are given in the title of each plot. Impact parameters with positive displacements lie on 
the red-shifted side of the rotational velocity gradient. The dashed lines show the predicted line profiles in the absence of rotation for 
comparison. 
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Figure 16. Simulated HCO+(J = 4 — > 3) and CS(J = 5 — > 4) centroid velocity contour plots (in kms ^) produced from the canonical 
infall model with no rotation. The centroid velocity was calculated over the velocity range —2.0 to +2.0 kms~^. 
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Figure 17. Simulated HCO+(J = 4 — + 3) and CS{J = 5 — > 4) centroid velocity contour plots (in kms ^) produced from the canonical 
infall model with a solid body rotation of f2sini=30kms^-'^ pc^"*^, calculated over the velocity range —2.0 to +2.0 kms^-'^. 
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Figure 18. Centroid velocity measurements of the predicted HCO^"(J = 4 — > 3) (open circles) and CS(J = 5 — > 4) (crosses) line profiles 
for the canonical infall model with a solid body rotation of f!sini=30kms~^ pc~^. The left hand panel shows a least-squares straight 
line fit to the centroid velocities measured along a line perpendicular to the rotation axis (i.e. in the direction of maximum rotational 
velocity gradient) passing through the centre of the cloud, while the right panel shows the centroid velocities along the rotation axis. 



© 0000 RAS, MNRAS 000, 000-000 



Modelling protostellar infall YI 



NGC 1333-IRAS2 



40 






20 




V(km/s) 

Figure 19. Spectra observed towards NGC1333-IRAS2. The ver- 
tical dashed line is at a velocity of 7.7 km s^^. An extra lOK base- 
line shift has been added to each successive spectrum. 



Date 



Reference 



1995 February 9-10 Feb95 

1995 July 26-28 Jul95 

1996 March 17-22 Mar96 

1997 January 30 Jan97 



Table 2. JCMT observations and their references. 

CS centroid velocities lie behind the overall velocity gradient 
on both sides of this gradient. 

This is due to the larger assumed FWHM beam size in 
the calculation of the CS(J = 5^4) profile (19.3" com- 
pared to 13.6"), which tends to increase the contribution of 
the brighter inner regions (which have lower rotational ve- 
locities) to the line profile. The centroid velocity along the 
axis parallel to the rotation axis in the right hand panel 
clearly shows the 'blue bulge' signature discussed above, al- 
though in this case the magnitude of the associated centroid 
velocity shift is rather small. 

The formal least-squares straight line fit to the cen- 
troid velocity gradient perpendicular to the rotation axis is 
23.5 ± 6.8kms~^ pc~^, and the intercept is — 0.018 km s~^. 
The actual angular velocity used in the model was Q. sin i = 
30kms~^ pc~^, so the actual and measured values agree to 
within la. The slightly low measurement for the centroid 
velocity gradient is mainly due to the lagging behind of the 
CS centroid velocities mentioned above, as can be seen from 
the figure. Nevertheless, this example shows that the cen- 
troid velocity gradient gives a reasonable estimate of the 



Source 


R.A.(1950) 


Dec. (1950) 


D 


Name 


h m s 


o / // 


(pc) 


N1333-IRAS2 


03 25 49.9 


+31 04 16 


220 


Serpens SMMl 


18 27 17.3 


-1-01 13 23 


300 


Serpens SMM4 


18 27 24.7 


+01 11 10 


300 


Serpens SMM3 


18 27 27.3 


+01 11 55 


300 


Serpens SMM2 


18 27 28.0 


+01 10 45 


300 



Table 3. Co-ordinates of the objects observed and their assumed 
distances. 



actual projected velocity gradient across a cloud, at least 
for the case of solid body rotation. 



5 OBSERVATIONS 

The observations were carried out at the James Clerk 
Maxwell Telescope (JCMT), using receivers RxB3i (Cun- 
ningham et al. 1992) and RxA2 (Davies et al. 1992), with 
the Digital Autocorrelation Spectrometer (DAS) backend. 
Details of the observing runs are given in Tables 2-4. Point- 
ing and focus checks were performed roughly once per hour 
during each run. The focus was found to be very stable 
throughout, and the pointing accuracy was ~ 2 arcsec. Ob- 
servations of standard sources were taken several times per 
night for each transition observed. We estimate the abso- 
lute calibration is accurate to within ~25%, although the 
relative calibration between different observations is much 
better than this. Throughout the observations the DAS was 
configured with the optimum frequency resolution of 95kHz 
per channel, giving a usable bandwidth of 125MHz (equiva- 
lent to < O.lkms"^ at the frequencies observed). 

Most of the pointed observations were made using po- 
sition switching. Frequency switching was used only for 
some of the HCO^ and H^^CO^ observations, and to check 
for emission-free reference positions during position switch- 
ing. Beam-switching was used only occasionally, for sources 
which were known to have very compact emission. Compar- 
ison of identical observations using each of these observing 
modes gives excellent agreement in each case - see Matthews 
(1996) for discussion of observing modes. Both grid mapping 
and raster mapping modes were used to make spectral line 
maps. In some of the maps, the noise is greater towards 
the edge of the map than at the centre, due to a reduced 
integration time per map point at these positions. In the 
individual raster maps, integration times of 4-6 seconds per 
cell were used, as recommended by the JCMT Users Guide 
(Matthews 1996). Longer integration times per point were 
obtained by repeating the raster. Data from the JCMT data 
archive were also used to supplement our observations where 
possible, as indicated in Table 4. 

Data reduction was carried out using the SPECX pack- 
age (Padman 1990). Baseline calibration was performed by 
subtracting a fit through sections of baseline placed on ei- 
ther side of the spectral feature of interest, carefully chosen 
to avoid suppression of low level wing emission. Frequency 
switching tends to produce curved or sinusoidal baselines. 
Curved baselines were removed by subtracting a second or- 
der polynomial fit in the vicinity of the line. This procedure 
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Molecule 


Transition 


Frequency 


FWHM 


Sv 


-^sys 


Reference 






(GHz) 


(") 


(kms-i) 


(K) 




HCO+ 


J = 3^2 


267.557625 


18.3 


~0.25 


~ 500 


archive 




J = 4^3 


356.734248 


13.6 


0.080 


~ 800 


Jul95/Mar96 


Hi3cO+ 


J = 3 ^2 


260.25548 


18.8 


~0.25 


~ 500 


archive 




J = 4-^3 


346.998540 


14.1 


0.082 


~ 800 


Jul95/Mar96 


CS 


J = 5-^4 


244.935606 


19.3 


~ 0.25 


~ 500 


Feb95/Jul95/Jan97 




J = 7-^6 


342.88294 


14.2 


0.083 


~ 800 


Feb95/Jul95/Mar96 


CO 


J = 3^2 


345.795979 


14.1 


0.082 


~ 800 


Mar96 


Ci80 


J = 2-^ 1 


219.56032 


22.1 


~ 0.25 


~ 500 


archive 


Ci80 


J = 3^2 


329.33050 


14.8 


0.087 


~ 1500 


Mar96 



Table 4. Table giving details of the submillimetre line transitions observed using the JCMT. The line frequencies were taken from the 
Lovas spectral line catalogue (Lovas 1992). The observing runs during which each transition was observed are indicated in the last column 
(c.f. Table 2). Also indicated are transitions which were obtained from the JCMT data archive. The CO(J = 3 — > 2) and CS(J = 7^6) 
transitions were observed simultaneously, by selecting a frontend frequency which placed one transition in each sideband. 



was checked on several occasions by making observations 
of the same object in both frequency-switched and position- 
switched mode, giving very good agreement between the hne 
profiles. 

The spectra are calibrated in terms of the main beam 
temperature scale, apart from the maps, where the corrected 
antenna temperature, TJ, is used. The integrated intensity 
and centroid velocity, position-velocity and channel maps 
were also produced using the standard SPECX map-making 
facilities. To improve the appearance of the contouring in 
these maps, the observed grids were first interpolated onto 
a grid with twice the spatial sampling rate. 

The relation between centroid velocity Vc over a given 
velocity interval in terms of Vi and TJ . , the velocity and an- 
tenna temperature of the i'th velocity channel, and the sum 
over the rich velocity channels lying in the specified velocity 
interval (Narayanan, Walker & Buckley, 1998) and axj^ the 
rms noise in a single channel, assumed to be constant for all 
the channels in the velocity window, is given by: 






+ 






2^^ E, "' 



(E.rij-(E,«^Tij 



This equation is used to derive the error bars for the centroid 
velocity measurements presented below. 



6 NGC1333-IRAS2 

6.1 Previous Observations 

NGC 1333 is a reflection nebula associated with the L1450 
dark cloud, which lies at a distance of ~220pc (Cernis 1990), 
in the Perseus molecular cloud complex (Sargent 1979). To 
the south of the reflection nebula lies a ~ 450 Mq molecular 
core, with a central density and temperature of 10'*cm~"^ and 
~ 18K respectively (Warin et al., 1996; Lada et al., 1974). 
This is a site of highly active low and intermediate mass star 
formation, evidenced by the large concentration of T-Tauri 
and Herbig Ae/Be stars (Lada, Alves & Lada 1996; Aspin, 
Sandell & Russell 1994), Herbig-Haro objects (Bally, Devine 
& Reipurth 1996), and bipolar jets and outflows (Hodapp & 
Ladd 1995; Liseau, Sandell & Knee 1988) in and around the 
core. Loren (1976) originally suggested that the star forma- 
tion in NGC1333 is being triggered by the collision of two 



dense molecular clouds. More recently, Warin et al. (1996) 
have argued that the morphology of the region supports a 
sequential star formation scenario, where outflows produced 
by one generation of stars compress the surrounding gas, 
which triggers collapse and further star formation. 

Jennings et al. (1987) listed nine compact IRAS sources 
associated with the NGC1333 cloud, flve of which drive 
bipolar molecular outflows (Liseau et al. 1988). Of these, 
IRAS 2 and IRAS 4 have no associated on-source optical 
or near-infrared emission. A full submillimetre continuum 
survey of NGC 1333-IRAS 2 was carried out by Sandell 
et al. (1994) using the JCMT. A strong compact peak 
was found at all wavelengths, with low-level emission at 
800pim extending northwest and southeast of the peak, in 
a 'flattened bar' morphology. Sandell et al. (1994) derived 
an envelope mass of ~ 0.8 M0 for the compact submil- 
limetre source. The total far infrared luminosity of IRAS2 
was estimated to be 17L0by Jennings et al. (1987) for a 
distance of 220pc. An upper limit to its total luminosity 
of 26L0was found by Ward- Thompson et al. (1996), who 
quoted the ratio of bolometric to submillimetre luminosity 
to be I/BOL/isuBMM < 130, below the Class threshold of 
200 (Andre et al., 1993). 

Two bipolar outflows have been associated with IRAS2, 
observed in CO ( J = 3 ^ 2) emission (Sandell et al. 1994), 
several transitions of CS (Langer, Castets & Lefloch 1996; 
Sandell et al. 1994), and 2.12^m shock-excited molecular 
hydrogen emission (Hodapp & Ladd 1995). The outflows are 
oriented nearly perpendicular to each other on the plane of 
the sky. However, since they typically only affect the extreme 
line wings, they do not alter the results of our modelling. 



6.2 New data 

Figure 19 presents spectra taken at the position of 
NGC1333-IRAS2 (hereafter IRAS2). Both the HCO+(J = 
4^3) and CS(J = 5^4) spectra show asymmetric dou- 
ble peaked line proflles, with the blue peak stronger. The 
absorption dip in the HCO^(J = 4 — > 3) line profile is 
particularly strong, and coincides in velocity with the peak 
of the rarer isotopomer lines. Also visible in the main line 
HCO^ and CS spectra are broad, low-level wings, which 
are probably tracing the outflowing gas associated with this 
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Figure 20. Grid map of T\ for NGC1333-IRAS2 in the HCO+{J = 4-^3) transition in units of K. 
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Figure 21. Grid map of T\ for NGC1333-IRAS2 in tfie CS(J : 



■ 4) transition in units of K. 



source. The peak of the C^*0(J = 2 — + 1) line is shifted 
towards the blue relative to the H^^CO"*'(J = 4^3) peak. 
Since H'^'^CO'''(J = 4 ^ 3) is expected to be the optically 
thinner of the two lines, we interpret this as showing that 
the C^®0 line is partially optically thick in this object. 

Several properties of the spectra agree with the pre- 
dicted qualitative signatures of infall in protostellar en- 
velopes (e.g. Myers et al. 1995; Zhou et al. 1993). For opti- 



cally thin lines, radiative transfer models predict a symmet- 
ric single peaked line profile, centred on the systemic veloc- 
ity. Such a profile is shown by the H^ CO^{J = 4^3) line, 
suggesting a systemic velocity of 7.7kms^^. The asymmet- 
ric double-peaked line profiles skewed towards blue velocities 
seen in both the HCO+(J = 4^3) and CS(J = 5^4) 
spectra are predicted by radiative transfer models of in- 
falling envelopes for optically thick lines. The velocity of the 
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self-absorption minimum may coincide with the systemic ve- 
locity, or be red-shifted with respect to it, depending on the 
optical depth of the line and the velocity and density struc- 
ture of the envelope. 

The C^*0(J = 2 -> 1) line profile, which is slightly 
asymmetric, and blue-shifted with respect to the systemic 
velocity, is also consistent with infall model predictions for 
a line with small but non-negligible optical depth. The only 
disagreements with the qualitative expectations of a pure 
infall model are the strong high velocity wings seen in several 
of the spectra, which we ascribe to the outflow (infall models 
do predict high velocity wings, but at a lower level than 
observed here), and the very small red-shift of the peak of 
the CS( J = 7^6) line profile, with respect to the inferred 
systemic velocity. 

Mardones et al. (1997) obtained on-source spectra of 
IRAS2 in CS(J = 2^1), N2H+(JfiF2 = loi ^ I02) 
and H2CO (Jk_iKi = 2i2 -^ In), using the IRAM 30-m 
telescope. The CS line is single-peaked and approximately 
symmetrical, peaking at the systemic velocity (7.7kms^^), 
and the N2H"*" line has an infall type asymmetrical profile. 
Within the infall scenario, this would suggest that the CS 
line is more optically thin than the N2H"'" line, contrary to 
expectations. The H2CO line profile observed by Mardones 
et al. (1997) has a complicated, multiple-peaked structure, 
with a very broad red-shifted wing, suggesting that this line 
is sensitive to the outflow. The line core is symmetrical and 
skewed to the red, which contrasts with the blue-skewed pro- 
files shown in Figure 20. The complete set of spectral line 
observations at the position of this object do not, there- 
fore, paint a wholly consistent picture. However, they are 
in reasonable agreement with the qualitative expections for 
an infalling envelope, despite the fact that outflow emission 
is probably contributing significantly to the line centres in 
many of the observed transitions. 

The spatial dependence of the observed line profiles 
may be able to shed light on this question. Grid maps 
of IRAS2 in HCO+(J = 4^3) and CS(J = 5-^4) 
emission are presented in Figures 20 and 21 respectively. 
In both maps the emission is centrally peaked, and the 
strength of the self-absorption feature decreases away from 
the peak. The off-source spectra in the maps show an in- 
teresting variety of profiles. In the 1100"*" map there is a 
clear trend towards broader linewidths from south to north, 
which is also discernable in the CS map. The CS map spec- 
tra show broader and relatively stronger wing emission than 
the HCO"*" spectra, and may therefore be more sensitive 
to outflowing gas. In a detailed submillimetre spectral line 
study of the NGC1333-IRAS4 proto-binary source, Blake 
et al. (1995) came to a similar conclusion, finding that the 
CS abundance is considerably enhanced in the outflow, and 
HCO^ appears to be less depleted than neutral species in 
the dense protostellar envelope. 

North of the map centre, a dramatic reversal of the line 
asymmetry in the HCO"*' map is seen, and the CS spectra 
show weaker blue-asymmetries. One possible explanation for 
this is that the emission from the northern (red-shifted) out- 
flow lobe is enhancing the redshifted emission, perhaps due 
to an interaction with ambient gas at this position. 

Centroid velocity maps calculated from the HCO^ and 
CS grid maps are shown in Figure 22. The interpretation 
of centroid velocity maps requires some care, since the cen- 



troid velocity may be quite sensitive to the velocity range 
over which it is calculated. To illustrate this, we show 
HCO+(J = 4^3) and CS(J = 5^4) centroid veloc- 
ity plots calculated for two different velocity windows, one 
covering the entire line profile, including the high velocity 
wings (3.7-11.7kms~^), and one covering mainly the line 
centre (6.2-9.2 kms~^). All of the plots show a clear veloc- 
ity gradient along an axis with a position angle between 0° 
and 20°. The velocity gradient is in the same sense as the 
north-south bipolar outflow suggesting a possible physical 
connection with it. The east-west outflow may be respon- 
sible for the small east-west velocity gradient apparent in 
the centroid velocity plots calculated over the wider veloc- 
ity window. 

The 'wide window' centroid velocity plots clearly show 
a much stronger velocity gradient than the 'narrow window' 
plots, particularly in the CS(J = 5^4) maps. This is 
consistent with the idea that this gradient is caused by the 
north-south bipolar outflow, and also supports the sugges- 
tion that CS is a more sensitive outflow tracer than HCO^ . 
The agreement in the direction of the centroid velocity gradi- 
ent for both the narrow and wide velocity windows suggests 
that the outflow even affects the line profile shapes in the 
line centres. This adds a serious complication to the inter- 
pretation of the line centre profiles in terms of infall models, 
since it is then extremely difficult to separate reliably the 
contribution of the envelope and outflow to the line profile. 

We now examine a possible alternative explanation for 
the velocity gradient seen in the centroid velocity maps, and 
the reversal of line asymmetry seen in the HCO^(J = 4 ^ 
3) map — rotation of the protostellar envelope. Rotation 
will tend to produce a centroid velocity gradient across the 
envelope, and radiative transfer models which include both 
rotation and gas infall predict reversals of the infall asym- 
metry at certain off-centre positions (e.g. Zhou 1995), al- 
though the degree of reversal in the HCO^ spectrum north 
of the map centre in Figure 21 is larger than typically pre- 
dicted by such models. Since the line wings of our CS and 
HCO^ observations are clearly affected by the outflow, we 
concentrate on centroid velocities calculated over the line 
centre only. 

Figure 23 shows linear fits to the combined CS(J = 
5 -> 4) and HCO+(J = 4^3) centroid veloci- 
ties calculated over the velocity window 6.7-8.7 km s^^, 
along north-south and east-west axes through the source. 
The fit to the north-south centroid velocity gradient is 
7.7(±0.9) X lO^^kms"^ 



arcsec 
this is equivalent to 7.2±0.8kms^^""^^ 



^ For a distance of 220pc, 
pc^^. The corresponding 



numbers for the east- west velocity gradient arc — 8.0(±0.8)x 
lO^^kms^^arcsec"^ and 7.5(±0.7) x 10"^ kms^^pc"^ The 
interpretation of centroid velocities is complicated by optical 
depth effects, since both the transitions are optically thick in 
the line centre. Furthermore, infall motions in the envelope 
will tend to skew the profiles of optically thick lines towards 
bluer velocities. 

One indication of this is the fact that the centroid 
velocities along the east-west axis are blue-shifted by ~ 
0.1 km s^^, on average, with respect to the systemic veloc- 
ity. The measured centroid velocity gradient for an optically 
thick line is therefore an unreliable measure of the actual ve- 
locity gradient in the object, and a proper treatment requires 
full 3-D radiative transfer modelling, or the use of optically 
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Figure 22. Comparison of centroid velocity maps of NGC1333-IRAS2 calculated over tlie velocity range 3.7— 11.71{rns^^(left) and 6.2- 
9.2kms~^(riglit), to cover tlie wliole line profiles and the line centres respectively. The top and bottom panels show HCO+(J = 4^3) 
and CS(J = 5 — > 4) centroid velocity plots respectively. The greyscaling and contour spacing (0.05 km s~^) are identical over all of the 
maps. The first light contour in each case marks a centroid velocity of 7.7 km s^^. 



thin lines to eliminate optical depth effects. However, for 
the purposes of the present discussion we assume that our 
centroid velocity gradient measurements give a reasonable 
indication of the actual line-of-sight velocity gradient across 
the object due to rotation. 

Using this assumption a lower limit for the binding mass 
can be estimated by hypothesising that the rotation is cen- 
trifugally balanced by gravity, and that the rotation axis lies 
in the plane of the sky. This gives a lower limit, because the 
actual rotation rate may be greater if the rotation axis is 
inclined out of the plane of the sky, and magnetic fields and 
pressure gradients may also contribute to the support of the 
envelope. Furthermore, if part of the envelope is infalling, 
as some of the observations suggest, then at some level the 
gravitational forces must dominate over the centrifugal sup- 
port. 

Assuming spherical symmetry, the mass M interior to 
a radius r is then given by: 



M = 



rv(rf r'^^{rf 



(31) 



G cos i G cos i 
where r is the distance from the centre of rotation, v{r) is 



the line of sight velocity at r (relative to the systemic ve- 
locity), f2(r) is the angular velocity of the rotation at r, G 
is the gravitational constant, and i is the inclination angle 
of the rotation axis from the plane of the sky. Using Equa- 
tion 31 with i = 0, we find a lower limit of 0.4 Mofor the 
enclosed mass inside a 30 arcsec (6600AU) radius. Sandell et 
al. (1994) estimated a total dust and gas mass in the IRAS2 
protostellar envelope of 0.79 M0, and to this must be added 
the mass of the central protostar. 



This analysis therefore does not contradict the hypoth- 
esis that the line-centre centroid velocity gradient is tracing 
rotation about an east-west axis projected onto the plane of 
the sky. If this were correct, the most likely model would be 
that the east-west bipolar outflow is driven by the Class 
source, and the protostellar envelope is rotating about this 
outflow axis. The older north-south outflow would then have 
to be driven by a separate, more evolved, and as yet unob- 
served object. This seems unlikely, and we consider it more 
likely that the velocity gradient is tracing the north-south 
outflow rather than rotation. 
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Figure 23. Plot showing weigtited least-squares linear fits to the centroid velocities along north-south (left) and east-west (right) 
axes through the position of the NGC1333— IRAS2 submillimetre continuum peak (at the map origin). The circles and crosses denote 
HCO+(J = 4^3) and CS(J = 5^4) measurements respectively. The centroid velocity was calculated over the line centre (6.7— 
8.7kms-l). 
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Figure 24. HCO+(J = 4^3) and R^'-^CO+{J = 4^3) (bold) spectra of the Serpens sources SMMl-4. 



7 SERPENS SMMl-4 

7.1 Previous observations 

The Serpens molecular cloud is a low- to intermediate-mass 
star-forming region which has been the subject of intensive 
observational study (e.g. Eiroa 1991; Eiroa & Casali 1992; 
White, Casali & Eiroa 1995, hereafter WCE; Hogerheijde et 
al. 1998; Testi & Sargent 1998; Testi et al. 2000; Williams 
& Myers 2000; Davis et al. 2000). The distance to the re- 



gion has been the source of some debate, with estimates 
ranging from 250 parsecs (Chavarria et al. 1988), to 750 
parsecs (Zhang et al. 1988), based on spectroscopic paral- 
laxes of stars believed to be physically related to the cloud. 
More recent studies (de Lara et al. 1991; Straizys, Cernis & 
Bartasiiite, 1996) have identified mis-classifications of the 
spectral types of some of the stars used in the previous 
distance determinations, and suggest an actual distance of 
300 ± 30 parsecs. 
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Eiroa & Casali (1992) carried out a comprehensive near- 
infrared survey of the region and positively identified 51 
young stellar objects associated with the cloud. They es- 
timated a lower limit of ~ 450 pc~^ for the average stellar 
density over a region of ~0.5 parsecs diameter, and a star 
formation efficiency in the range 8-28 per cent, suggesting 
that a bound cluster may be forming. The peak of the star 
formation activity is centred close to the Serpens reflection 
nebula (SRN), which is a diffuse optical and near-infrared 
nebula (Gomez de Castro, Eiroa & Lenzen 1988; Eiroa & 
Casali 1992) powered by the luminous pre-main sequence 
star SVS2 (Strom, Vrba & Strom 1976). 

Maps in ammonia (Ungerechts & Giisten 1984), C^'^O 
(WCE), CS (McMulIin et al. 1994b) and far-infrared emis- 
sion (Zhang et al. 1988; Hurt & Barsony 1996) show the 
gas and dust in the central region to be concentrated in 
two cores, aligned approximately northwest-southeast. Es- 
timates of the gas density and temperature in the cores 
lie in the range 10* — lO^cm"^ and 20-30K respectively, 
with higher temperatures and densities towards some of the 
embedded sources. The ~7OM0 northwestern core contains 
the well studied far-infrared source Serpens FlRSl (Har- 
vey et al. 1984), which is associated with a highly colli- 
mated, high velocity (~ 300 kms~^) radio continuum jet 
(Curiel et al. 1993; 1996), and a linear near-infrared fea- 
ture, approximately aligned with the northwestern lobe of 
the radio jet (Eiroa & Casali 1989). This source was labelled 
SMMl (Casah, Eiroa & Duncan 1993, hereafter CED), and 
is the strongest of the submillimetre continuum sources in 
the cloud. 

The more massive southeastern core includes the reflec- 
tion nebula, and contains a small cluster of compact submil- 
limetre continuum sources (CED), two of which (SMM2 and 
SMM4) have no near-infrared counterparts. Barsony (1997) 
used a maximum correlation algorithm to enhance the spa- 
tial resolution of the IRAS maps of this region, and derived 
upper limits for the far-infrared luminosities of the submil- 
limetre sources. SMMl-4, and S68N (McMullin et al. 1994b) 
were identified as Class sources. Hurt, Barsony & Woot- 
ten (1996) have observed each of these sources in the high 
density molecular tracer H2CO. In the 'R2CO{Jk_iKi = 
3o3 -^ 202) transition, four out of the five sources showed 
asymmetrical line profiles suggestive of infall. The remaining 
source, SMMl, showed a symmetrical double-peaked profile. 
Gregersen et al. (1997) included Serpens SMMl-4 in their 
HCO^ survey of Class sources, and found strong infall 
signatures towards SMM2 and SMM4. 

The CO (J = 2 ^ 1) outflow map of WCE 
shows possible bipolar outflow associations with SMMl 
(southeast-northwest) and SMM4 (approximately north- 
south), and perhaps SMM2 (southeast-northwest). Estab- 
lishing the presence or absence of outflows driven by the 
other Serpens objects is hindered by the considerable con- 
fusion from overlapping outflow lobes. 



7.2 New data 

Figure 24 shows the HCO+(J = 4^3) and \i^^CO+{J = 
4^3) spectra observed towards the Serpens Class sources 
SMM 1-4. Of the four sources, SMM2 and SMM4 show hue 
profiles most suggestive of infall. SMM4 in particular shows 
the classic signature of a double-peaked self-absorbed main 



line profile skewed heavily towards the blue, with an absorp- 
tion minimum redshifted relative to the peak of the optically 
thin isotopic line. The SMMl main-line profile shows a deep 
absorption minimum, close to the peak velocity of the iso- 
topic line, and a very mild blue asymmetry. Strong high 
velocity wings are also visible, which are probably tracing 
the energetic bipolar jet and outflow driven by this source. 
The HCO+(J = 4^3) spectrum observed towards SMM3 
does not show any sign of self-absorption, but a strong blue- 
shifted wing, probably tracing an outflow, is clearly seen. In 
the following we concentrate primarily on the sources SMM2 
and SMM4, as the best protostellar infall candidates in the 
Serpens cloud. 

Figure 25(a) shows the on-source spectra observed to- 
wards Serpens SMM2. Blue asymmetric line profiles are seen 
in most of the transitions, although CS(J = 7^6) emission 
was not detected above the ~0.2K noise of the observation. 
There is a remarkable similarity between line profiles of the 
main line HCO^ transitions and their isotopic counterparts, 
although the signal to noise ratio in the H^''CO"'"( J = 4-^3) 
line profile may be too low to attach much weight to this. 

Figure 26(a) shows a grid map of Serpens SMM2 in 
the HCO'''(J = 4 — > 3) transition. This map was obtained 
from the JCMT data archive, and the observing runs dur- 
ing which these data were taken are described in WCE. The 
HCO'^(J = 4^3) emission is seen to peak northwest of the 
continuum peak. However, the bright line profiles observed 
towards the north-western corner of the map include a con- 
tribution from the outer envelope of SMM4, which lies at an 
offset of (—25,-1-50), and shows extended bright emission in 
HCO^(J = 4^3)- see Figure 27. This may in part explain 
the similar velocity structure of the isotopic and main line 
profiles discussed above, if the emission arises from overlap- 
ping clumps with different velocities along the line of sight. 
There is no reliable way to separate the overlapping contri- 
butions to the line profiles, and we therefore do not discuss 
this source any further. 

The spectra observed towards the position of the Ser- 
pens SMM4 submillimetre continuum peak are shown in Fig- 
ure 25(b). The optically thin lines - H"CO+(J = 3-^2) 
and H"CO+(J = 4 -> 3) - peak at a velocity of 7.95 kms^^ 
The C^*0(J = 3^2) line peaks at 7.7kms'\ suggesting 
that there is a non-negligible optical depth in this line, as we 
have found previously. There is good qualitative agreement 
between the line profiles shown and typical radiative transfer 
model predictions for an infalling envelope (see the discus- 
sion above of the NGC1333-IRAS2 on-source spectra). 

Both the main line HCC^ transitions show the asym- 
metric double-peaked line profile skewed towards blue veloc- 
ities. We also note that the CO(J = 3^2) line also shows 
this type of profile. However, CO line profiles are not as 
robust infall indicators as HCO-f, since they suffer greater 
confusion from foreground material and outfiowing gas. The 
isotopic HCC^ line profiles shown are also in qualitative 
agreement with model predictions, being single-peaked, and 
peaking between the velocity of the absorption minimum 
and the blue-shifted peak of the main-line profile. 

Further qualitative agreement between model predic- 
tions and observations is shown by three isotopic CO lines. 
Models predict that as the optical depth of the observed 
transition increases, the line profiles progress from single 
peaked and symmetric, to single-peaked and skewed to the 
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Figure 25. (a) Spectra of various transitions observed tovifards the position of the Serpens SMM2 continuum peak. The dashed line 
indicates a velocity of 7.9 km s^-*^. (b) Spectra observed towards the Serpens SMM4 continuum peak. The dashed vertical line indicates 
a velocity of 7.95 km s^^. 
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Figure 26. (a) Grid map of Serpens SMM2 in the HCO+(J = 4^3) transition (taken from the JCMT data archive). The origi 
the position of the submillimetre continuum peak, (b) Grid map of Serpens SMM4 in the HCO+(J = 4 — » 3) transition, centred 
submillimetre continuum peak. 



in is at 
on the 



blue, to double-peaked with the blue peak brighter than 
the red (Myers et al., 1995). This is just the progression 
observed in the CO line profiles as the isotopic abundance 
increases. The CS line profiles represent the intermediate 
optical depth stage, being single peaked and skewed to the 
blue. Taken as a whole, therefore, the set of line profiles 



shown in Figure 25(b) are remarkably consistent with the 
qualitative expections for an infalling envelope. 

Figure 26(b) shows a grid map of the SMM4 enve- 
lope in the HCO^(J = 4^3) transition. The strongest 
emission is found at the origin and at the grid position 
immediately north of the origin. Blue-skewed line profiles 
are evident in most of the spectra shown, although double 
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Figure 27. (a) Comparison of spectra in tlie northwest (—14, +14) and southeast (+14,-7) corner of the HCO+(J = 4 — » 3) grid map 
of Serpens SMM4 shown in Figure 26(b). (b) Ccntroid velocity map derived from the HCO^"(J = 4 — » 3) grid map of Serpens SMM4, 
over the velocity range 5.0-ll.Okms^^. The contour spacing is 0.05 kms^^. 



peaked line profiles are only found in the central and west- 
ern map columns. Comparison of the spectra in the south- 
east and northwest corners of the map reveals a shift in the 
peak velocity of about Ikms^^towards the blue along the 
southeast- northwest direction - see Figure 27(a). A centroid 
velocity map calculated from this grid map is shown in Fig- 
ure 27(b), from which a southeast-northwest centroid veloc- 
ity gradient of ~ 30kms~^pc~^ is measured. Also apparent 
in this plot is a 'blue-bulge' encroaching onto the red-shifted 
side of the velocity gradient from the blue-shifted side. As 
discussed above, this is a predicted signature of infall in the 
presence of rotation. 

A larger area HCO+(J = 4^3) raster map of SMM4 
was also made, and is shown in the left hand plot of Figure 28 
as an integrated intensity map. The HCO^(J = 4 — > 3) 
emission is seen to be strongly enhanced in the region of 
SMM4 and is well resolved in our beam (the transition from 
dark to light contours in the figure marks the half-maximum 
level). The emission is approximately centred around the 
offset position (0,+5) on the map. Towards this point the 
integrated intensity fiattens off into a plateau, and there is 
a shallow minimum at the central point itself, surrounded 
by horseshoe-shaped peak. The lower contours reveal an ex- 
tension of the peak towards the northeast, and a flattened 
linear feature oriented northwest-southeast. 

The right hand panel of Figure 28 shows the centroid 
velocity map calculated from the raster map. The southeast- 
northwest velocity gradient is again very apparent. Compar- 
ison of integrated intensity and centroid velocity maps shows 
that the region of brightest integrated HCO^(J = 4^3) 
emission coincides with the blue-shifted part of the velocity 
gradient, which makes the symmetric appearance of the plot 
about the line rest velocity all the more remarkable. The axis 
of the steepest velocity gradient (P.A.~ —40° ± 5°) is close 
to (but not exactly coincident with) the axis of the elongated 
feature in the integrated intensity map (P.A.~ —50° ± 3°). 

The steepest velocity gradient is —4.1 (±0.15) x 



10^^ kms^^arcsec^^, which implies a gradient of 28 ± 
1 km s^^pc^^ , for a distance of 300pc. The centroid velocities 
along this axis clearly shift towards bluer velocities at posi- 
tions near to the central peak. This 'blue-bulge' signature, 
parallel to the rotation axis of a rotating, infalling proto- 
stellar envelope, was first predicted by the radiative transfer 
analysis of Walker et al. (1994). 

As before, we estimate the lower limit on the mass as- 
suming the observed velocity gradient is due to rotation and 
the rotation axis lies in the plane of the sky. We consider the 
mass enclosed within the limits of the steep velocity gradi- 
ent between position offsets of —20 and +10 arcsec from 
the HCC^ peak, and obtain a lower limit of ~ 1.9 M© (for 
a distance of 300 parsecs). Hurt and Barsony (1996) esti- 
mated a gas mass of 3 M0 surrounding SMM4, derived from 
dust continuum emission, which, along with the mass of any 
embedded protostars, is sufficient to provide the centripetal 
force required to bind the rotation, as long as the rotation 
axis is close to the plane of the sky. Nevertheless, if our 
analysis is correct, centrifugal effects are likely to have a 
very significant influence on the equilibrium and dynamics 
of the envelope gas. The elongation seen in the integrated 
HCO^(J = 4 — > 3) map, which is aligned approximately 
perpendicular to the inferred rotation axis, may well in fact 
be caused by centrifugal flattening of the outer envelope (e.g. 
Pudritz & Norman 1986). 

The alternative possible explanation for the observed 
velocity gradient in the HCO"^(J = 4^3) maps could be 
outflowing gas. Gregersen et al. (1997) detected a velocity 
gradient in the same sense and along the same axis in their 
HCO+(J = 3 -» 2) map of SMM4, but interpreted it in 
terms of outflow rather than rotation. SMM4 has previously 
been associated with a north-south bipolar outflow (WCE), 
but not aligned with the gradient we observe. Examination 
of the channel maps in our own data likewise reveals no 
evidence for an outflow in this direction. Thus we believe 
our rotation hypothesis to be the more likely explanation. 
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Figure 28. HCO"'"(J = 4 — > 3) integrated intensity (left) and centroid velocity (right) maps of Serpens SMM4, over the velocity interval 
4-12krns~^. The peak contour of the integrated intensity map is 18.6Kkms^^, and the contour spacing is 0.6Kkms~'^(~ la). In the 
centroid velocity map, the minimum and maximum contours levels are at 6.8 and 8.8kms~^respectively, and the contour spacing is 
O.lkms"^. The rest velocity of the source (~ 7.95 km s~^) lies in the transition interval between the dark and light contours. The cross 
indicates the effective FWHM beam size. 
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Figure 29. Model fits (dashed lines) to the on-sourcc spectra (solid lines) observed towards NGC1333-IRAS2. The assumed systemic 
velocity is 7.75 kms^-*^, which has been subtracted from the velocities in the observed spectra (see text for details). 



8 MODELLING THE DATA 
8.1 NGC1333-IRAS2 

Figure 29 shows our best model fit to the HCO+ and CS 
line profiles observed towards NGC1333-IRAS2. We ap- 
proached the modelling by first searching for a fit to the 
HCO"'"(J = 4^3) profile, choosing realistic physical pa- 
rameters. We aimed to achieve consistency in our model fits 



with previously determined properties of the objects, such 
as the envelope mass and ambient gas temperature. Once 
a suitable fit was found, the same model was used to find 
the best fit to the 11^^00+ (J = 4^3) profile, but in this 
case the only free parameter was the [H^^CO+]/[H^^CO+] 
ratio, which we required to lie in the range 50-100. Depend- 
ing on the quality of this fit, the model was either rejected, 
or tested further against the CS line profiles. 
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Parameter 



Value 



^;.(r<6000AU) [-0.29 (g^jji^) ^ 

I n on / r \ ^ 



0.29 ( 



N 0.15 



6000 AU ) 



Vr (r > 6000 AU) km s" 
nn^CrOOOAU) 2.46 x 106[0.35 (g^^)" 
^0.65(^^^-°-''i- 



kmf 



\ 900 AU / 

nH2(r>900AU) 2.46x106(3^5^2:^) 
T(r<2400AU) 17 (^^gl^) -"■' - 



)-"-^^] 



T{r > 2400 AU) 
A?;tb(FWHM) 

'''^HCO + 
12C/13C 



17 K 
0.60kms-i 
6.0 X 10-9 
6.0 X 10^9 
90 



Table 5. The parameters used in the model of NGC1333-IRAS2. 



We attached less weight to the CS profiles in the model 
fitting, since CS tends to be more affected by outflow emis- 
sion than HCO+ (e.g. Blake et al. 1995). Langer et al. (1996) 
mapped the entire NGC1333 core region in the CS(J = 5 ^ 
4) , (J = 3 ^ 2) and (J = 2 ^ 1) transitions using the 
IRAM telescope, and found that the CS emission extended 
over a very wide area. The CS profiles may therefore be sig- 
nificantly affected by emission and absorption from the gas 
in the foreground core. 

The parameters used in the model fit are listed in Table 
5. The outer radius of the cloud was set at 10 000 AU, and an 
innermost shell radius of 15 AU was used. A distance to the 
cloud of 220 pc was assumed (Cernis 1990). The total enve- 
lope mass within the outer radius of the cloud inferred from 
the molecular hydrogen density profile is 3.7 M0. The values 
of Tinf used in the infall velocity and density profiles were 
6000 AU and 900 AU respectively, and Ocft = 0.29 km s~^was 
used in both cases. No acceptable fit could be found which 
used identical values of rinf for both profiles. In other words, 
we could not find a solution where the density and veloc- 
ity profiles taken together are consistent with the singular 
isothermal sphere model. 

The fit to the HCO+(J = 4-^3) line profile is very 
good between —1 and -|~1 km s^^ (relative to the systemic 
velocity of 7.75 kms~^). At larger velocities in both the red 
and blue wings of the line, the model significantly underesti- 
mates the emission. The most likely origin for the excess high 
velocity emission is from one or both of the bipolar outfiows 
associated with this source - the clear break in the slope 
of the observed line profile at ~ ±1.3 km s"^ also suggests a 
separate origin for the high velocity emission. We have not 
tried to estimate how much the outflow contributes to the 
line profile at lower velocities, so this is another source of 
uncertainty in the model. 

The H^^CO+(J = 4^3) model fit agrees weU with the 
peak velocity and intensity in the observed line, although 
the model profile has slightly overestimated the line width. 
The width of the H"C0+(J = 4^3) line in the model is 
mostly determined by the systematic and turbulent veloci- 
ties near the centre of the cloud, since the warm dense gas 
in this region has the strongest intrinsic emission, and for 
an optically thin tracer this emission is unattenuated by the 



outer envelope on its way to the observer. The asymmetry 
in the main line profile is mainly sensitive to the systematic 
and turbulent velocities at larger radii. If the infall velocity 
at large radii is chosen to produce agreement with the asym- 
metry in the HCO^(J = 4 — > 3) profile, then the inferred 
velocity at small radii always produces too much broadening 
in the H^^CO+( J = 4^3) line, as a result of the v ex r"^/^ 
velocity law. 

The largest discrepancy between the observations and 
the model is in the CS(J = 5^4) line, where the model 
predicts a substantially higher intensity in the line core than 
is observed. Better fits to the overall intensity of the CS line 
(but not the line shape) could be found by reducing the CS 
relative abundance, although this then underestimated the 
intensity of the CS(J = 7 — > 6) line. This might be explained 
if the CS relative abundance may vary with radius, but we 
did not pursue this. The value of the CS relative abundance 
used in the model fits was simply optimised to give the best 
agreement with the CS(J = 7^6) line. 

Figure 30 shows an alternative fit to the HCO^(J — 
4^3) and H"CO+(J = 4^3) line profiles obtained by 
setting the infall velocity to zero inside a radius of 800 AU, 
keeping all other parameters the same. The main line profile 
is affected remarkably little by this dramatic alteration to 
the inner velocity field (although the fit is slightly worse), 
which underlines the point that each observed transition 
constrains only a limited domain of parameter space in the 
model envelope. The CS line profiles (in which the CS rela- 
tive abundance has been lowered by a factor of 4 relative to 
the original model fit) are still not well fitted by this model. 
In both Figures 30 & 31 the predicted H"C0+(J = 4^3) 
and CS(J = 7^6) line profiles are almost identical to each 
other, and we could not produce simultaneous fits to the 
velocity widths of both observed profiles. Since CS emission 
shows a greater tendency for contamination by the outflow, 
we are inclined to attribute the larger velocity width on the 
CS(J = 7^6) line to the outflow, and therefore we tend 
to prefer the former fit described above. 

In Figure 31, we compare the off-centre IICO^(J — 
4-^3) line profiles with the model predictions. Because 
of the assumption of spherical symmetry, the predicted line 
profiles depend only on the impact parameter of the line of 
sight relative to the origin, and not on its direction. The 
decrease in line intensity with increasing impact parameter 
is reasonably well matched by the model. However, the line 
shapes at the off-centre positions show less good agreement. 
With the exception of the spectrum at the (0,-1-15) position, 
the off-source line profiles are all skewed towards the blue. 
The model profiles are much more symmetric, as a result of 
the fall-off in infall velocity at large radius, which is again a 
consequence of the v ex r^^" infall velocity law. It is very 
difficult to accommodate the strong high velocity emission 
seen in several of the off-source positions - most prominently 
at (0,-1-15) and (-1-15,-15) - within a plausible infall model, 
particularly given the highly asymmetric distribution of the 
high velocity emission in both position and velocity, which 
cannot be fitted by a spherically symmetric model. 

We explored the possibility above that the centroid ve- 
locity gradient measured across the NGC1333-IRAS2 enve- 
lope using HCO+(J = 4^3) and CS(J = 5^4) spectral 
line maps could be explained by rotation, and concluded 
that although rotation could not be ruled out, outfiow was 
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Figure 30. Model fits using the same model parameters as in Figure 30, but with the infall velocity set to zero inside a radius of 800 AU, 
and a factor of 4 reduction in the CS relative abundance. 



the more likely explanation. On the basis of the results in 
Figure 31, we believe that the north-south bipolar outflow, 
and not rotation, is responsible for the centroid velocity gra- 
dient. 

The mass of the central protostars can be approximately 
estimated from the infall velocity profile used in the model 
fit, applying the freefall equation M* = rv'^ /2G to the inner 
part of the profile. Using this relation we obtain an estimate 
of ~ O.3M0. However, this must be regarded as somewhat 
uncertain. The density and velocity profiles adopted in this 
model are inconsistent with the SIS model, despite individ- 
ually having the SIS model forms. The infall radius used in 
the velocity profile is a factor of ~ 7 greater than the infall 
radius used in the density profile, indicating that the infall 
is more well developed than a SIS model envelope with a 
similar density profile. 

Hydrodynamical (e.g. Foster & Chevalier 1993) and 
magneto-hydrodynamical (e.g. Ciolek & Mouschovias 1994) 
simulations of the evolution of dense cloud cores during the 
pre-stellar phase predict that by the time the density profile 
approximates closely to a singular distribution, substantial 
infall velocities have already developed. Notwithstanding the 
uncertainties in our model fit, our results support this pic- 
ture of a non-static initial condition for collapse. 



8.2 Serpens SMM4 

Figure 32 shows the model flts to the CS and HCO"'" on- 
source spectra observed towards Serpens SMM4. The ap- 
proach to the modelling was similar to that used for 
NGC1333-IRAS2, starting with the HCO+(J = 4^3) 
line, and then making adjustments to produce the optimal 
agreement with the remaining lines. To produce a fit to the 
IICO^(J = 4-^3) line required the infall radius in the SIS 



model form of the velocity proflle to be set well beyond the 
outer radius of the cloud, i.e. the entire protostellar envelope 
must be inf ailing. 

If the infall radius is set inside the outer radius of the 
cloud, then absorption from the static gas beyond this radius 
tends to produce an absorption dip close to the systemic 
velocity, in disagreement with the observed proflle. As the 
SIS model velocity proflle is inappropriate in this case, we 
used a simple r^^ " power law, which is the asymptotic form 
of the inner part of the SIS model proflle, to model the 
velocity profile instead. 

The parameters used in the model are given in Ta- 
ble 6. The outer radius of the model envelope was taken as 
10 000 AU, and we used an innermost shell radius of 62 AU. 
The assumed distance to the source was 300 pc. The total en- 
velope mass implied by the model is ~5.1 M©, roughly con- 
sistent with the estimate of ~3 Mq derived from the IRAS 
fiuxes by Hurt & Barsony (1996). 

As with NGC1333-IRAS2, the infall velocities in the 
model fit are much greater than would be predicted by the 
SIS model given the density proflle. From the infall veloc- 
ity proflle we estimate the mass of the central protostellar 
system to be M* = r«^/2G — 2.1Mq. The temperature pro- 
file is normalised at a considerably higher level than in the 
IRAS2 fit, which would suggest that SMM4 is considerably 
more luminous. 

A high accretion luminosity for this object is also sug- 
gested by the large central protostellar mass and mass ac- 
cretion rate (~ 6 x lO"'^ M© yr"^). Hurt & Barsony (1996) 
estimated the luminosity of SMM4 as 9 Lq , after deconvolv- 
ing the HIRES IRAS images of the Serpens region to obtain 
the far-infrared fiuxes. Jennings et al. (1987) estimated the 
luminosity of IRAS2 as 17Lq, also using IRAS fiuxes, so 
there are clearly large uncertainties on these estimates. The 
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Figure 31. Comparison of the oflf-centre HCO+(J = 4^3) line profiles observed towards NGC1333-IRAS2 (solid lines) with the model 
predictions (dashed lines). The offsets (in arcsee) from the IRAS2 position arc given in the title of each plot. 



Parameter 



Value 



-l-95(Too^Au)"*kms-i^ 



2000 AU ; 
-0.64, 



\ cm 
cin^' 



Vr(r) 

nH2(r- < 2000AU) 1.48 x 10*5 [0.35 ( 

+0-65(200^; 

nH2 (r > 2000 AU) 1 .48 x lO^ ( 505^ ) 

T(r< 3800 AU) 2^^^,r_)-'-'-' k 

T(r > 3800 AU) 25 K 

At)tb(FWHM) O.SOkms^i 

Xhco+ 3.5 X 10-a 

Xcs 3.5 X 10-^ 

i^C/i^C 90 



Table 6. The parameters used in the model fit for Serpens SMM4. 

parameters used in our niodel fits are inconsistent witli tliese 
luminosities, and for our model to be correct, then either the 
SMM4 luminosity calculated from the IRAS data, or our 
assumed protostellar radius, would have to be a significant 
underestimate. 

The predicted H^'^CO^ lincwidths are considerably 
broader than the observations, as we found with NGC1333- 



IRAS2, although the peak intensities are in very good agree- 
ment. Again, we suggest this may be indicating deviations 



of the infall velocity profile below the 



.-1/2 



free-fall pro- 



file in the centre of the cloud. We found above that the 
SMM4 envelope shows strong evidence for rotation (f2 sin i ~ 
28 kms~^pc~^), and suggested that the flattening seen in 
the spatial distribution of integrated HCO^(J = 4^3) 
emission may be caused by centrifugal support of the outer 
envelope. Thus centrifugal braking may retard the infall ve- 
locity near the centre of the cloud. 

The calculated CS(J = 5^4) and CS(J = 7^6) 
line profiles match the observed line intensities reasonably 
well, but the peak velocities of the calculated profiles are less 
red-shifted than the observations. The enhanced blue-shifted 
wings seen in both of the observed CS spectra strongly sug- 
gest that outflow contributes to the emission in these lines. 
As discussed above, CS is believed to be enhanced relative 
to HCO+ in bipolar outflows (Blake et al. 1995). 

Figure 33 shows the flts to the off-centre HCO^(J = 
4-^3) line profiles towards SMM4 from our best-fit model. 
We included a solid body rotation of 28kms^^ pc^'^ in the 
model (assuming rotation axis in plane of sky), consistent 
with the above centroid velocity gradient across SMM4. As 
the HCO+(J = 4^3) spectra at the (0,0) and (0,+7) 
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Figure 32. Model fits to the CS and HCO"'" on-source spectra observed towards Serpens SMM4. 



positions have approximately equal intensities, we assumed 
that the centre of the cloud lies at the offset (0,+3.5). Al- 
though the model profiles somewhat underestimate the ob- 
served lines, we are encouraged by the overall agreement of 
the trends in the observed lineshapes across the envelope. 
It is likely that a better fit to the line intensities requires a 
shallower density profile than predicted by the SIS model. 

In contrast to IRAS2 we find little evidence for out- 
flow emission in the off-centre line profiles, and essentially 
all of the high velocity emission can be accounted for by the 
infall. We consider these results to provide convincing evi- 
dence that both infall and rotation are taking place in the 
the Serpens SMM4 envelope. The magnitude of the rota- 
tional velocity is large enough for centrifugal effects to have 
a significant infiuence on the gas dynamics. 



9 SUMMARY 

In this paper we have discussed the physical constraints 
which can be used to limit the number and range of pa- 
rameters needed to define the spherically symmetric models 
in the Stenholm radiative transfer code. We have studied 
the dependence of the predicted line profiles on a number 
of model parameters. Increasing the infall velocities in the 
envelope tends to produce stronger asymmetry between the 
blue and red-shifted peaks of double-peaked line profiles, 
without strongly afi'ecting the separation between the peaks. 
The value of the turbulent velocity dispersion is the most 
significant factor in determining the velocity separation be- 
tween the peaks. We carried out a preliminary investigation 
into the efl'ect of adopting different power law dependences 
of the turbulent velocity width on the gas density and found 
the predicted line profiles to be extremely sensitive to the 
value of the power law exponent. 



We have studied the effect on the predicted line pro- 
files of super-imposing a solid-body rotation onto the infall 
velocity field. The effect of the solid-body rotation on the on- 
source double-peaked line profiles was to 'smooth out' the 
double-peaked structure. The centroid velocities and the line 
shapes of the off-centre profiles were found to be significantly 
altered by the solid-body rotation. On the red-shifted side 
of the rotational velocity gradient, a reversal of the asymme- 
try in the line profiles was seen. We measured the centroid 
velocity gradient predicted by our model for a solid-body ro- 
tating envelope, which was found to slightly underestimate 
the actual rotational velocity gradient used in the model. 

We selected the objects NGC1333-IRAS2 and Ser- 
pens SMM4 for detailed modelling using the Stenholm ra- 
diative transfer code, confining our search to the SIS model 
forms for the density and velocity profiles in the modelling, 
although we did not require these profiles to use the same 
SIS model parameters in the models. Our best model fits to 
both the NGC1333-IRAS2 and Serpens SMM4 data used 
velocity profiles which were more strongly developed than 
the SIS model would predict, given the fit to the density 
profile. This is consistent with theoretical predictions that 
at the instant a central protostar is formed, the gas in the en- 
velope has already accelerated to substantial velocities (e.g. 
Foster & Chevalier 1993; JCH), which contrasts with the 
static initial state assumed by the SIS model, and is one of 
the principal distinguishing features between the SIS model 
and other collapse models. 

Our fits to the H^^CO^ line profile observations of both 
objects slightly over-estimated the profile widths. We could 
produce better fits to these lines by lowering the infall veloc- 
ity in the centre of the cloud, although this affected the fits to 
the main line profile slightly. The infall velocity may there- 
fore be retarded from the r~^ " free-fall relation towards the 
centre of the cloud, possibly by centrifugal braking. In the 



© 0000 RAS, MNRAS 000, 000-000 



ModellinQ vroto stellar infall 31 



(7,14) 



(0,14) 



(-7,14) 











(0,7) 





( 


A ' 


5 


/ 


\ ( 





f^' 


/"Vt^ 



(-7,7) 



V(km/s) 




(7,0) 



(0,0) 







(-7,0) 







1 / \ 


- 


5 


- // \ 


- 





^''' n 


%m 



c 

V(km/s) 



(7.-7) 



(0,-7) 





-4-2024 -4-20 

V(km/s) V(km/s) 




-2 

V(km/s) 



Figure 33. Comparison of tlie oflf-centre HCO^(J = 4^3) line profiles observed towards Serpens SMM4 (solid lines) with the model 
predictions (dashed lines). The offsets (in arcsec) from the IRAS2 position arc given in the title of each plot. In the model spectra we 
have included a solid body rotation of 28 km s~^ pc^^ , about an axis passing through the centre of the cloud at a position angle of +40° . 
The centre of the cloud was assumed to lie at the offset (0,+3.5). 



case of IRAS2, the discrepancies between the model predic- 
tions and the observations for the ofF-source profiles could 
not possibly be accommodated within a plausible model of 
infall and rotation, and we attribute them to the outflow. 

The model fits to the off-centre HCO+(J = 4^3) line 
profiles of Serpens SMM4 included the solid-body rotation 



inferred from the HCO"'"(J = 4 — > 3) centroid velocity gra- 
dient in the calculation. The predicted and observed profiles 
showed good agreement in the overall line shape although 
the model slightly underestimated the line intensities. We 
thus conclude that Serpens SMM4 is an infalling and rotat- 
ing protostar. 
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